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SECTION I 


INTRODUCTION 


This final report has been prepared by the General Electric Company, Aerospace Electronic 
Systems Department, Utica, New York under contract NAS8-28516. The report documents the 
results of a thermal spreading resistance data generation technique study. The method developed 
is discussed in detail, illustrative examples given, and the resulting computer program is in- 
cluded. 
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SECTION n 
BACKGROUND 


A. GENERAL 

"Thermal spreading resistance" is defined as the conductive thermal resistance between a 
source region and a sink region in a solid where the geometry is such as to preclude one dimen- 
sional heat flow. 

Knowledge of thermal spreading resistance is needed in two aerospace engineering areas. 
These are the thermal design of electronic components or equipments and in the prediction and 
control of thermal contact resistance. 

1. Importance To The Design Of Electronic Components and Equipments 

The thermal analysis of a power semiconductor or integrated circuit can be reduced to 
the problem of determining the appropriate spreading and bonding thermal resistances. As an 
example, the problem of calculating the junction-to-case thermal resistance of a semiconductor 
bonded to a substrate which is bonded in a metal case will be considered. Figure 1 illustrates 
this problem. 


JUNCTION 



Figure 1. Semiconductor in an Integrated Circuit 

Heat is generated in a region of known size, the junction region of the semiconductor. 
The first, and most significant, spreading resistance of interest occurs between the junction and 
the opposite face of the silicon chip. The next thermal resistance of interest is that across the 
bond between chip and substrate. It is of significance that these thermal resistances are not in- 
dependent although many thermal designers, under the pressures of a design schedule, have 
treated them as such. The thermal conductance of the bond proper can vary several thousand- 
fold depending on the use of a metallic or nonmetallie bonding material. The resistance to heat 
flow between the semiconductor chip bond region and the rear of the substrate represents a 
second spreading resistance, etc. In a typical integrated circuit package the entire bottom 
region of the substrate would not be available as a sink for a single semiconductor chip due to the 
presence of other heat dissipating chips. It is usually possible to estimate the effective sink re- 
gion on the rear of the substrate from considerations of symmetry or because it exceeds dimen- 
sions which appreciably affect the thermal spreading resistance. In those few cases where inter- 
actions must be considered, the key analytical tool is superposition; Green’s function approach 
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may also be employed to advantage. For example, see reference 1 and the discussion beginning 
on page 37 of this report. 

The importance of being able to predict thermal spreading resistances in single and 
multi- layered material in the evaluation of the thermal design of semiconductor or integrated 
circuits has been shown. Spreading thermal resistances are important in other electrical devices 
such as phased array antenna elements, Peltier coolers, Seebeck generators and many devices 
which utilize conductive heat transfer. 

B. PREDICTION AND CONTROL OF THERMAL CONTACT RESISTANCE 

The resistance to heat flow between two mating (touching, as in a joint) pieces of metal 
is called thermal contact resistance. When the actual microscopic regions of contact between 
two mating surfaces are examined, it is found that metal-to-metal contact occurs in small dis- 
crete regions where the asperities or microscopic protuberances make contact. References 2 
and 3 describe this model of contact in great detail. Figure 2 illustrates this contact model. 



Figure 2. Microscopic View of the Joint of Contacting Pieces of Metal 


The heat flow to and from a region of asperitic contact into the contacting proper 
is seen to be of the "spreading" type. In fact, the effective thermal contact resistance of any 
contact may be considered as the sum of the parallel microscopic spreading resistances in the 
contacts themselves. References 2 and 3 above deal largely with isentropic contacts in which the 
thermal conductivity within the bodies of both contacts is uniform. 

Analysis has shown that the bulk of the spreading resistance occurs close to the region 
of actual asperitic contact and that the spreading resistance in any region varies inversely with 
the thermal conductivity of the material. Figures 3 and 4 illustrate the first of these points. 
Figure 3, drawn to scale, shows the equipotential lines about a circular contact region each 
drawn to show one-tenth of the total spreading resistance between the circular source region and 
the body of a very large contact. It is seen that half of this total resistance occurs within one 
contact radius from the circular contact or source region and 80 percent occurs within three 
contact radii. Figure 4 illustrates these relationships. Figures 3 and 4 are taken from 
reference 4. 

The thermal conductivity of contact close to the surface is of such importance that 
even a thin 45 Angstrom thick layer of oxide on an aluminum contact can contribute measurably 
to the thermal contact resistance of an aluminum contact. This has been shown by Gale, 
reference 4. 

Mikic and Carnasciali, reference 5, have utilized the above principle to enhance ther- 
mal contact conductance by plating materials of higher conductivity on the contacting faces of a 


r 
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Figure 3. Temperature Profiles Described by Holm's Equation for Isothermal Circular Source 

on a Semi-infinite Slab 
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DISTANCE FROM SURFACE -CONTACT RADII 

Figure 4. Percent of Total Constriction Resistance for a Single Isothermal Circular Source on 
a Semi- infinite Slab as a Function of Distance into Body of Contact 


metallic joint. They have attempted an analysis of spreading resistance from a circular contact 
into a contact composed of two layers of materials with different conductivities. An exact 
boundary value solution of this basic problem has proven too difficult as no mathematical function 
has been found which will satisfy the boundary conditions between the plating and the body 
materials. 

Professor C. J. Moore, Jr. * in his discussion printed at the end of reference 5 felt 
this two layered spreading resistance problem could best be handled by a ’Veil-conditioned 
finite difference computer code. " Mikic and Carnasciali then question the economic feasibility 
of such calculations. 

v/-\. 

Attempts by the author of this study to solve the two layer thermal spreading resistance 
problem using a finite difference approach utilizing Gauss- Seidel iteration have shown the cost of 
digital computer calculation to be great. 
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SECTION in 


A. GENERAL 


THEORY 


The governing differential equation for the thermal spreading resistance problem is 
Poisson’s equation. For those spreading resistance problems that are two-dimensional^ or may 
be reduced to two-dimensional problems, the equation is: 





(1) 


Consider a rectangular field subdivided into rectangular subregions as illustrated in 
Figure 5. The heat balance equation describing the heat flow among element m, n and its four 
principal neighbors is : 


(T 

m,n 


-T 


m,n+l' 


H 


m,n 


+ (T - T , ) 

m,n m-l,n 


V +(T -T J 
m,n m, n m , n- 1 


H 


m,n- 


1 


where : 


+ (T -T )V 

m,n m+l,n' m+l,n 


*m,n 


( 2 ) 


T is temperature 

q’” heat generated per unit volume 
x,y,z are spatial coordinates 

H,V are horizontal and vertical conductances, respectively 

q heat generated in mode m,n 
m , n 


The convention for the horizontal and vertical conductances used is shown in Figure 6. 


Each of the following observations below will be helpful in understanding the discussion 
which follows: 


(1) When any temperature t m is known (e.g. , as a boundary condition), it will 
affect equation m,n by yielding a term qm } ii, which is subtracted from the 
right hand side of equation (2) where q m n ’’ is: 


^m,n 


’ = T (H + V + H i + V , ) 

m,n m,n m,n m,n-l m+1, n 


(3) 


^Associate Professor of Mechanical and Aerospace Engineering, North Carolina State University, 
Raleigh, N.C. 

^The method developed is applicable to three-dimensional problems as will be shown later in the 
report. 
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(2) If the original field is divided into M rows and N columns, and further if 

M = N, then: 

(a) There will be linear equations. 

(b) There will be not more than unknowns (fewer if some temperatures 
are initially prescribed). 

(c) There can be as many different and distinct nodal conductances as there 
are interconnections between nodes. 

Now, if the system of linear finite difference equations is written in matrix form (taking the 
nodes of Figure 5 into consideration) from left to right, top row to bottom row, as in reading 
English, a coefficient matrix results that has a pattern characteristic for field problems described 
by Poisson s or LaPlace’s equations. This pattern is illustrated in Figure 7. 

ft was noted by Karlqvist (Reference 6) that the matrix in Figure 7 may be partitioned as 
shown. It can be seen that each of the submatrices is N x N and the coefficient matrix is N 2 x N 2 
where the original finite element matrix was N x N in size 


B. DERIVATION OF AN EFFICIENT 
OF EQUATIONS 


TECHNIQUE FOR EXACT SOLUTION OF THIS SYSTEM 


Defining the sub- matrices shown in Figure 7 as follows: 


B] C 1 0 0 0 


V 


V 

&2 &2 C 2 0 0 


T 2 


q 2 

0 A3 83 c 3 0 


t 3 

zz 

°3 

° ° A 4 B 4 C 4 


T 4 


°4 

0 0 0 a 5 b 5 


_ T 5_ 


_ Q s_ 


Figure 8. System Of Submatrices In Matrix Notation 

Expanding the partitioned matrices (Figure 8) into a system of equations, 
having normalized each equation with respect to the diagonal element: 


T 1 

VVss 

0 

0 

0 = 

Y 1 ^ 

VVl 

T 2 

VVs 

0 
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V 1 ^ 

0 

vV 2 

T 3 

V 1C 3 T 4 
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0 

0 

4 4 3 
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B 4 _1 Q4 

0 

0 

0 

VS T 4 

T 5 = 
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Upon redefining constants in the following manner: 

^2 “ "® 2 * B 2 ^2 anc * B 2 ^2 ~ ^ 2 * e * c ' 
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Figure 7. System of Equations in Matrix Notation 



the general equation has the form : 


-B.T. , + T. - A.T. ... = C. 
li-l 1 1 l+l l 

The first equation can be solved for T^: 

T i ■ C 1 + A 1 T 2 
and the ith for T. : 

l 

T. = C. + A.T. i + B.T. , 
l l li+l li-l 


( 4 ) 


(5) 


( 6 ) 


The goal is to find a recursion relationship built upon successive substitutions, which pro- 
vides a solution for the ith unknown in terms of the (i+l)th. That is: 


T. = A.’ T „ + B * 
l li+l i 


(7) 


Examining equation (5) for above, it can be seen that: 


A 1 = A x and 


The equation for T is: 
£ 


T 2 = C 2 +A 2 T 3 + B 2 T 1 


( 8 ) 


which, when written in terms of the equation for T^, becomes: 


T 2 ' 


I - B^ 


-1 


A 2 T 3 + 


i-b 2 a 


i r i 


C 2 + B 2 B l' 


The general coefficients found in this manner become 

1-1 


A i' = 


^Vi-I 


A. 


(9) 


and 


B. f = 

I - B.A. / 

-1 

B.B. / 

+ C. 

l 

i l-l J 


l l-l 

l 


Therefore 


T. = A. f T. « + B.' 
l l l+l l 


( 10 ) 


The temperature matrices (columns) are found starting at the Nth row by making 


T N = V 


( 11 ) 


A^ = 0 as a boundary condition results in modification of above [ see equation (3) ] . 


The system of equations has been solved by operating on 3\/N~ - 2 submatrices, each of 
which is the square root of the size of the original NxN coefficient. inversions of these 
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submatrices are required. The total number of multiplications (an indication of the effort) re- 
quired during solution is: 

No. of Multiplications = 3N 2 + N ?//2 - N + ' (12) 

This may be compared against other direct methods (see Ref. 7): 

Number of Multiplications 

Method Required during Solution 


Gaussian Elimination 

Jordan 

Doolittle 


1 xt 3 , tvt 2 1 ir 

3 N +N ' 3 N 

1 TU 3 XT 2 1 XT 

-hr + n - - N 


Cholesky ^N 3 +|n 2 +-jN 

Cornock's method (Ref. 8), a triangulation type, also makes use of the characteristic 
pattern of submatrices which results during a finite difference solution for fields described by 
Poisson’s equation. When the field properties are homogeneous and isentropic, Cornock's 
method is very powerful since only one of the above submatrices of order V N need be inverted. 
However, for the general solution of the nonhomogeneous field, the number of multiplications 
required is 


i! N 2 - — N 3/2 
2 W 2 N 


2N - 5 


(13) 


A serious drawback to Cornock’s method is that it does not lend itself to ready general program- 
ming for matrices of variable size as does the method described in this report. 

That equation (12) is indicative of the computer effort required for solution has been sub- 
stantiated in practice. Figure 9 shows the variation in cost realized in the solution of very large 
matrices using the method developed in this report. Also, a strong feature of this method is that 
very large systems of equations, e. g. , 2500, can be conveniently handled in a direct solution. 

The FORTRAN Y program contained in the Appendix was used on a Honeywell 630 computer 
in generating the dollar costs shown in Figure 9. Out-of-core storage of submatrices was utilized 
for very large systems. 


C. APPLICABILITY OF TECHNIQUE TO THREE-DIMENSIONAL PROBLEMS 

The technique discussed above is suited to the solution of field problems having three 
or more dimensions. Figure 10 illustrates the characteristic pattern of the coeffici ent m atrix 
for a three dimensional finite element array. It is seen that the submatrices are\3/ N in size. 
The same block tridiagonal pattern of submatrices is seen to occur as in the two-dimensional 
case so the derivation above for the technique of solution for two dimensional matrices is still 
applicable. Thus, although the BASIC and FORTRAN Y program presented later in this report 
are written for two-dimensional problems, little revision of thes e programs is required to handle 
three dimensional programs. Since the submatrices are \y N rather than the X^/lTin size, the 
technique is even more powerful for three dimensional problems. The number of multiplications 
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Figure 9. Cost of Computation vs Size of Coefficient Matrix 
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required for solution of the three-dimensional problem is a function of N as opposed to 
the two dimensional array where N is the order of the coefficient matrix. 
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SECTION rv 


ILLUSTRATIVE PROBLEM AND PROGRAM 


A sample thermal spreading resistance problem will be solved to illustrate the tech- 
nique presented in Section II. The computer program used in the problem is written in BASIC 
language. A general version of the same program written in FORTRAN Y is included in the 
Appendix. 

A. DESCRIPTION OF SAMPLE PROBLEM 

The thermal spreading resistance problem to be considered is depicted in Figure 11. Heat 
is uniformly generated in a plane circular region of radius a and flows to a circular sink of radius 
b both concentric with, and parallel to, the source region a distance H away in a conductive 
medium. The conductive medium is divided into two regions of different conductivity- An exact 
closed form solution of this problem has not been found. 


T MAX 



Figure 11. Mathematical Model for Spreading Resistance Nomographs 
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Ratios of geometries and conductivities used in the generation of sample data are sum 
marized below. All data is generated and presented in nondimensional parameters. 


Parameter 

Number of Values 

Parametric Values 

a/b 

6 

0.111, 0.222, 0.166, 0.551, 0.388, 0.5 

H/b 

6 

0.1, 0.2, 0.5, 1, 2, 5 

C/H 

5 

0.1, 0.2, 0.3, 0.5, 

K1/K2 

5 

0.01, 0.1, 0.2, 0.5, 1 


The above three-dimensional problem can be viewed as two-dimensional since all heat flow 
within the cylinder is in the axial and radial directions. Further, there is symmetry about the 
axis of the cylinder- 

Figure 12 shows the arrangement of finite elements used in the illustrative data generation 
program. 

B. CALCULATION OF NODAL CONDUCTANCES 

The convention used for the nodal conductances was that the vertical conductance V(m,n) 
associated with each node was that in the upward direction and the horizontal conductance ^ 
was that connecting with the node on the right. This is illustrated in Figure 13. 


(m,n+l) 


Figure 13. Convention for Nodal Conductances 

The calculation of conductances is straightforward for all modes except for those nodes of 
horizontal conductance lying on the axis of the cylinder, i.e., n = 1. The horizontal conductance 
of this inner node was approximated by using the exact solution of Jacob (Ref. 9) for two- 
dimensional heat flow within a cylinder having uniformly distributed internal heat generation for 
the difference between the mean temperature of the cylinder and the outside surface with radial 
flow. 




The conductances of all other nodes could be calculated in an exact manner using calculus. 
General expressions for the M, Nth nodal conductances in terms of M, N were developed and used 
in the sample problem to facilitate changing the program to allow the use of different numbers of 
finite elements. 
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C . DETAILED DESCRIPTION OF ILLUSTRATIVE PROGRAM 

Figure 14 shows the computer program for the illustrative program. The major steps in the 
program are described below . 


Line No. Description 

10-60 Dimensioning symbols will be defined as used in this program. 

The same symbol may be redefined several times. 

70 P = 7T 


80 Cl is that part of the horizontal thermal resistance of nodes 10, 1 

and 1, 1 between the nodal point and the surface of these nodes. 
(See discussion above concerning Jacob’s formula.) 

90 Expressions for entire horizontal resistance between nodes 10, 1 

and 10, 2 as well as 1, 1 and 1, 2. 


110 


H is the horizontal conductance. 


160 

200 

210 

290, 320 

380-460 

470 

480-250 

530-651 
680-690 
720-801 
810-971 
990, 1000 
1010-1220 
1010 
1080 


V is the vertical conductance to node above. 

Generalized expression for horizontal conductance for most 
modes, m, n. 

Generalized expression for vertical conductance for most modes, 
m, n. 

Entering adiabatic boundary conditions at top and sides . 

A9 is the radius of the heat generating region. 

Calculation and print of a/b (see Figure 11). 

Entering uniform heat input in the circular region described by 
A9. Use is made of the area dependence share by vertical con- 
ductance and heat input. 

Setting of H/b (see Figure 11). 

Adjustment of horizontal and vertical conductance for H/b. 

Setting of C/H (see Figure 11). 

Establishment and assignment of values for K1/K2 (see Figure 11). 
Optional printout o£ H (m> r) and V (m> n) 

Generation of coefficient matrix. 

M is row number of physical nodal pattern. 

Sets subdiagonal and superdiagonal in the coefficient matrix (line 
arrays W immediately either side of the main diagonal) to zero 
(see Figure 15d) . 

(Continued on page 25) 
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RECOMM GE RECOMM GE RECOMM GE RECOMM 


1 PRINT'* A/9 "l"H/8 **J M C/H **I M KI/K2 

2 PRINT 

10 DIM H< 12# I 2) # VC 1 2# t 2) # T< 12# 1 2) # QC 10# 10) 

20 DIM XC10* I0)*YC10# 10># WC 10# !0)*2C10# !0> 

30 DIM AC 10* 10) *8( 10# 10>* DC 10# 10>#EC 10* 10) 

40 DIM FC 10# 10>* G< 10* 10)# !< 10# 10>*J< 10* !0)*KC 10# 10) 
SO DIM LC10* 10) *M( 10* 10>*NC 10* 1O>*0< 10* 10)* P< 10* 10) 
60 DIM R< 10* 10)*SC 1 0* 10)* U( 10# 10) 

70 P* 3» 14159265 
90 Cl». 125/P 
90 R1»C!*(L0GC2))/C2*P> 

100 R2-R1/2 

101 F0R N4-4 T0 5 

102 FOR N3* 1 T05 

103 F0R N2-1 T0 6 

104 FOR N 1 « 1 T09 

110 HC 1 > 1 >*HC 10# 1 )» 1/R1 
120 F0R M»2 T0 9 
1 30 HCM* l )«2*HC1# 1 ) 

140 NEXT M 
15OP0R H* 1 T0 10 
1 60 VCM* l )«P/2 
170 NEXT M 
1 90 F0R M»! T0 10 
1 90 F0R N-2 T0 10 

200H(M#N> * 1 /( l/<49P)*L0GCN/CN-t )) ) 

210 VCM#N)-P*4*CN» 1 ) 

220 IF M* 1 • 1 THEN 240 
230HC l#N)*.5*HC l#N) 

240 MOO#N)».5*H( 10»N) 

250 VCM# 10>»P*C4*N-5>/2 

260 NEXT N 

270 NEXT M 

290 FOR M«1 T0 10 

290 HCM* 10) *0 

300 NEXT M 

310 FOR N* 1 T0 10 

320 VC 1 *N > *0 

330 NEXT N 

380 IF N 1 * 1 THEN 2660 

390 IF N | *5 THEN 2660 

400 IF Nl=7 THEN 2660 

410 IF N1 *9 THEN 2660 

420 49® 1 

430 IF N I <2 THEN 470 
440 A9«2*Nl-l 
450 IF N 1 < 1 0 THEN 470 
460 49-18 

470 PRINT USING 471*49/18# 

471! **•## 

480 CC1#t)-VC2* 1)92.00000 

490 IF N 1 ® 1 THEN 530 

491 F0R N®2 T0 10 

492 CCN# 1 >«0 

Figure 14. Computer Program for Illustrative Example 
(Sheet 1 of 6) 
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COMM GE RECOMM GE RECOMM GE RECOMM GE 


493 NEXT N 

500 F0R N*1 T0 Nt 

510 CCN# |>-VK2jN>*2-0000 

520 MEXT « 

S30 HI * . 1 

540 IF N2<2 THEN 650 
550 H ! ■ • 2 

560 IF M2<3 THEN 650 
570 HI = • 5 

580 IF N2< 4 THEN 650 
590 H l * I 

600 IF N2<5 THEN 650 
610 Hl-2 

615 IF N2< 6 THEN 650 
620 H 1 » 5 

630 IF N2<7 THEN 650 
640 H ! « I 0 

650 PRINT USING 651#H1# 

651 I #M##I#.## 

660 FOR H=1 T0 10 
670 F0R N-! T0 10 
680 V<M,N>*V(M#N)/HI 
690 HCM#N)=H<M*N)*H1 
700 NEXT N 

710 NEXT M 
720 Cl- I 

730 IF N3<2 THEN 800 
740 C 1-2 

750 IF N3<3 THEN 800 
760 Cl * 3 

770 IF N3<4 THEN 800 
780 Cl - 5 

785 IF N3<5 THEN 800 
790 Cl- 7 

800 PRINT USING 801*C1/!0* 

801 * IIIIHM# 

8 1 0 K2-1 

820 IF N4«2 THEN 910 
830 K2«2 

840 IF N4< 3 THEN 910 
850 K2-5 

860 IF N4< 4 THEN 910 
870 K2-10 

880 IF N4<5 THEN 910 
890 K2-100 

910 F0R N»C I ♦ ! T0 10 
920 F0R N«1 T0 10 
930 V<N,N)-V<M,N>9K2 
940 H(M>N)«H<N*N)*K2 
950 NEXT N 
960 NEXT M 

970 PRINT USING 972*1/K2* 

972t 

980 G0 T0 1010 

Figure 14, Computer Program for Illustrative Example 
(Sheet 2 of 6) 
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MM GE RECOMM GE RECOMM GE RECOMM GE RE 


990 MAT PRINT Hi 
tOOO MAT PRINT VI 
1010 FOR M«1 TO 10 
1020 MAT X*ZER 
1030 MAT r*ZER 
1040 MAT W«ZER 
1050 FOR NM T0 10 
1060 IF M* 1 THEN 1090 
1070 IF N»10 THEN 1090 
1080 W(N#N-l>«W<N#N+!)=0 
1090 IF N< 2 THEN 1 1 10 
1100 W<N#N-1)=-H(M#N-|> 

1110 IF N> 9THEN 1 ! 30 
1 120 WCN,N* 1 >*-H<M#N> 

I 130 XCN#N)*-V<M+1#N) 

1 140 Y<N*N)»-V(M,N) 

1 150 C3*0 

I 160 IF N*t THEN 1 t80 
1170 C3«VKN#N-1) 

1180 IF N* 1 0 THEN 1200 
1190 C3»C3*W<N*N* !) 

1200 W(N#N>«C3*X(N/N)*rCN,N> 
1210 W(N>N)=-W(N»N) 

1220 NEXT N 

1230 MAT Z* INV<W> 


1240 

MAT 

w* z*x 


1250 

MAT 

T*<- 1 >*W 


1260 

MAT 

W=T*U ) 


1270 

MAT 

x* z*r 


1280 

MAT 

T=ZER 


1290 

MAT 

T» C - 1 > *X 


1300 

MAT 

X« T4 C I ) 


1310 

IF Mo 1 THEN 

1 350 

1320MAT A-Z*C 


1330 

MAT 

C=ZER 


1340 

MAT 

B-WPC1 ) 


1350 

IF M<>2 THEN 

1 380 

1360 

MAT 

C*X*< 1 ) 


1370 

MAT 

o« w* < n 


1 380 

IF M<> 3 THEN 

1 410 

1390 

MAT 

E°X* C l ) 


1 400 

MAT 

F«W*( 1 ) 


1410 

IF M<> 4 THEN 

1 440 

1420 

MAT 

G»X*< 1 ) 


1430 

MAT 

I * W4 ( 1 ) 


1440 

IF M<» 5 THEN 

1 470 

1 450 

MAT 

J=X*< | ) 


1 460 

MAT 

K=W*< 1 ) 


1 470 

IF Mo 6 THEN 

1 500 

1 480 

MAT 

L*X* f 1 > 


1 490 

MAT 

Matf*< 1 ) 


1500 

IF Mo 7 THEN 

1530 

1510 

MAT 

N«X* < 1 ) 


1520 

MAT 

0*W*C t > 



Figure 14. Computer Program for Illustrative Example 
(Sheet 3 of 6) 
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GE RECOMM GE RECOMM GE RECOMM GE RECO 


1530 IF M<>8THE N 1560 
1540 MAT P*X* < 1 > 

I 550 HAT O^W+C 1 > 

1560 IF H<*9TWEN 1590 
1570 MAT R»X*< 1 ) 

1580 MAT S=W+< 1 ) 

1590 IF M«10 THEN 1630 
1600 MAT T«ZER 
610 MAT T*X* Cl) 

620 MAT U*W*< ! ) 

630 NEXT M 

640 FOR M»1 T0 10 

650 FOR N»i T0 10 

660 W<M,N>*XCM,N>»Y(M,N>aZCM,N)®0 

670 NEXT N 

680 NEXT M 

690 MAT X-C*B 

700 MAT Y= I ON 

710 MAT Z*Y-X 

720MAT X*INV fZ> 


1 730 

MAT 

Z»X»D 

1 740 

MAT 

0«Z*C1 > 

1750 

MAT 

Z*C*A 

1760 

MAT 

c«x*z 

1770 

MAT 

X* E*D 

1780 

MAT 

Z=Y-X 

1790 

MAT 

XalNV(Z) 

1800 

MAT 

Z»X*F 

1810 

MAT 

F«Z*C1> 

1820 

MAT 

Z«E*C 

1830 

MAT 

E»X*Z 

1840 

MAT 

X« G*F 

1850 

MAT 

Z* Y-X 

I860 

MAT 

X* I NV( Z > 

1870 

MAT 

Z« X*I 

1880 

MAT 

I»Z*<1> 

1890 

MAT 

Z*G*E 

1900 

MAT 

G»X*Z 

1910 

MAT 

XaJ*I 

1920 

MAT 

Z«Y-X 

1930 

MAT 

X* I NV( Z> 

1940 

MAT 

Z»X*K 

1950 

MAT 

K»Z*C 1 ) 

1960 

MAT 

Z«J*G 

1970 

MAT 

J»X#Z 

1980 

MAT 

X«L*K 

1990 

MAT 

Z*Y-X 

2000 

MAT 

X*INV(Z> 

201 0 

MAT 

Z«X*M 

2020 

MAT 

M«Z*( 1 > 

2030 

MAT 

Z«L* J 

2040 

MAT 

L«X*Z 

2050 

MAT 

XaN*M 

2060 

MAT 

Z»Y-X 


Figure 14. Computer Program for Illustrative Example 
(Sheet 4 of 6) 
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COMM GE RECOMM GE RECOMM GE RECOMM 


2070 

MAT 

X»IMV<Z> 

208 0 

MAT 

Z*X40 

2090 

MAT 

0* Z4 ( 1 > 

2100 

MAT 

Z«N*L 

2110 

MAT 

X+Z 

2120 

MAT 

X*P*0 

2130 

MAT 

Z«Y-X 

2140 

MAT 

X*INV<Z> 

2150 

MAT 

Z = X*Q 

2160 

MAT 

0«Z*C1> 

2170 

MAT 

Z«P*N 

2180 

MAT 

P*X*Z 

2190 

MAT 

X«K*9 

2200 

MAT 

Z* Y-X 

2210 

MAT 

X* IMVCZ> 

2220 

MAT 

z*x*s 

2230 

MAT 

s*z*c n 

2240 

MAT 

z* n * p 

2250 

MAT 

r*x*z 

2260 

MAT 

X* T* S 

227 0 

MAT 

Z«Y-S 

2280 

MAT 

X* I MV CZ> 

2290 

MAT 

Z*X*U 

2300 

MAT 

U«Z*< 1 ) 

2310 

MAT 

Z*T 

2320 

MAT 

T*ZER 

2330 

MAT 

X* ZER 

2340 

MAT 

X=S*T 

2350 

MAT 

Y»ZER 

2360 

MAT 

Y*X* R 

2370 

MAT 

R-Y * ( 1 > 

2380 

MAT 

X a Q*Y 

2390 

MAT 

Y«X + P 

2400 

MAT 

P* Y* C 1 > 

2410 

MAT 

X»i*Y 

2420MAT Y*X^N 

2430 

MAT 

N*Y»CI> 

2440 

MAT 

X*M* Y 

2450 

MAT 

Y«X*L 

2460 

MAT 

L*Y* ( 1 > 

2470 

MAT 

X*K*Y 

2480 

MAT 

Y*X*J 

2490 

MAT 

J* Y* ( 1 ) 

2500 

MAT 

X* I *Y 

2510 

MAT 

Y 3 X*G 

2520 

MAT 

G* Y*< 1 ) 

2530 

MAT 

X=E* Y 

2540 

MAT 

y»x*e 

2550 

MAT 

E- Y*f 1 ) 


2560MAT X*D *Y 


Figure 14. Computer Program for Illustrative Example 
(Sheet 5 of 6) 
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RECOMM GE RE 


2570 MAT Y«X*C 
2580 MAT C»Y*Cn 
2590 MAT X«B*Y 
2600 MAT Y*X+- A 
2610 MAT A«Y*(I ) 

2630 R7“ AC 1 # 1 ) /A9 
2640 PRINT USING 26S0-R7 
2650 If###.## 

2660 NEXT N1 
2665 PRINT 
2670 NEXT N2 
2675 PRINT 
2680 NEXT N3 
2685 PRINT 
2690 NEXT N4 
2695 PRINT 
2700 END 


Figure 14. Computer Program for Illustrative Example 
(Sheet 6 of 6) 
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N,P,R,T 


M,0,Q,$ 






























line No. Description 

1100 Assigns values to first subdiagonal of the coefficient matrix W 

(see Figure 15d) describing the mth row of the physical nodal 
pattern. 

1120 Assigns values to first subdiagonal of the coefficient matrix W as 

in line 1100. 


1130 

1140 

1200, 1210 

1230 

1240 

1270 

1280 

1320-1630 

1640-1680 


Assigns values to the X submatrix (see Figure 15c). 

Assigns values to the Y submatrix (see Figure 15b). 

Enters elements down the main diagonal of the W submatrices. 
First submatrix manipulation statement. 

Normalizes the submatrix to the right of the W matrix diagonal 
(see Figure 15a). 

Normalizes the submatrix to the left of the W matrix 
(see Figure 15a). 

Empties T matrix. 

Assigns each of the submatrices, subdiagonal and superdiagonal, 
mapped by Figure 15a after normalization. Note that matrices 
X and W are functions of M, the row of the physical nodal pattern. 

Clears matrices W, X, Y, Z, so they may be redefined below by 
entering only nonzero elements. 


1690-2250 


2320 

2330-2610 


These steps calculate the recurrence coefficients A^’ and B^' ac- 
cording to equation (10) of Section II. As these are calculated, 
the superdiagonal and subdiagonal matrices of Figure 15a are 
sequentially redefined to be these recurrence coefficients. 

This statement enters the boundary condition that the temper- 
atures along the bottom row of the physical matrix are zero. 

Using the recurrence relationship developed in 1690-2250, 
equation (10) is used to calculate the temperatures, one row 
(of the physical nodal model) at a time. These (column) matrices 
of temperatures are calculated in the following order (see Fig- 
ure 15a) and with the following nomenclature: T, R, P; N, L - 
J, G, E, C, A. 


2630 R7 is the temperature of node 1, 1 divided by the radius of the 

heat source A 9, defined in statements 420-460. 

D. SAMPLE SOLUTIONS 


Figure 16 presents the sample solutions of the illustrative problem. t^AX is the 
average temperature nodal point 1, 1, see Figure 12. It is obvious that, for these solutions, 
increasing the number of nodes in the physical model would have the result of increasing the 
temperature in this hottest element. This was done and the results are discussed below. 
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Figure 16. Sample Solution of Illustrative Problem 
(Sheet 1 of 10) 
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Figure 16. Sample Solution of Illustrative Problem 
(Sheet 2 of 10) 
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Figure 16. Sample Solution of Illustrative Problem 
(Sheet 3 of 10) 
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Figure 16. Sample Solution of Illustrative Problem 
(Sheet 4 of 10) 
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Figure 16. Sample Solution of Illustrative Problem 
(Sheet 5 of 10) 
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Figure 16. Sample Solution of Illustrative Problem 
(Sheet 6 of 10) 
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Figure 16. Sample Solution of Illustrative Problem 
(Sheet 7 of 10) 
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Figure 16. Sample Solution of Illustrative Problem 
(Sheet 8 of 10) 
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Figure 16. Sample Solution of Illustrative Problem 
(Sheet 9 of 10) 
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Figure 16. Sample Solution of Illustrative Problem 
(Sheet 10 of 10) 
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Effect of Increasing the Number of Finite Elements on the Solution of the Illustrative 
Problem 

An exact closed form solution exists for the cylindrical spreading resistance problem 
in a medium of uniform conductivity. Kennedy (Ref. 10) shows that, as a/b -*• 0, 

^MAX-y Kl _ 

q” a 

Figure 17 shows this trend for the finite difference model. With 1600 nodes, (tMAX' t s) Kl/q” was 
calculated to be 0. 9788. With 2500 nodes, the nondimensional resistance dropped to 0. 9766; this 
reduction is attributed to the inexact treatment of the horizontal resistance of the nodes in the 
inner column using the equation of Jacob (Ref. 9) as described earlier. 




a = I 



0 500 1000 1500 2000 2500 

NUMBER OF NODES IN PHYSICAL MODEL 

Figure 17. Maximum Nodal Temperature for Minimum a/b as a Function of Number of Nodes 


It would appear that a nodal model of 900 nodes would be a near optimum number for 
generating solutions to this particular problem using the finite element approach. 
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SECTION V 


THEORY OF SUPERPOSITION OF SOLUTIONS OF SPREADING THERMAL 

RESISTANCE PROBLEMS 


A . GENERAL 

The general steady state heat-conduction equation in Cartesian coordinates, known as the 
Poisson equation, is given by 


where 


Sft . j^t_ £t_ 

ax 2 a y 2 6z 2 


+ F ■ o 


(14) 


t = t(x, y, z) = temperature (° F) 

K = thermal conductivity (taken to be independent of temperature and position) 
(Btu/hr-ft-° F) 

qi.» _ q..» y> z j _ internal volumetric heat source (Btu/hr-ft 3 ) 
x, y, z = Cartesian coordinates 

The generalized boundary conditions vary. For example, specified temperature: 


t(x, y, z) 


x b ,y b ,z b = 


(15a) 


where: 


f = specified function, 

x^, y b , z^ = values of x, y, z on the boundary (b) 
or, convection to a fluid: 


where 


bt 

Sn 


V y b’ % 


_h 

k 


t(x, y, z) 


on the boundary 



n = outward directed vector normal to the boundary 

2 

h = Newtonian convective film coefficient (Btu/hr-ft -°F) 
tj = bulk temperature of convecting fluid (°F) 


(15b) 


r - r- 
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It is convenient to rewrite equations (14), (15a), and (15b) using a temperature difference 
for the dependent variable that contains a reference temperature. This reference temperature is 
typically taken to be a specified boundary temperature, as in equation (15a), or the fluid temper- 
ature as in equation (15b). 


Hence, we define this temperature difference by 
U ~ * ^reference 

Equations (14), (15a), and (15b) can then be rewritten 

\2 --.2 \2 

du, du , d u , ^ 

~T + TT~ + 7T +G - 0 

dx d y 0 Z 

U{X ’Mx y ^ “ s ( Wb> 


du 
d n 


where 


x b> y b> *b 
q m /K 


= ~K 


(16) 


(17) 

(18a) 

(18b) 


These equations are linear as can be seen by observing that they contain no products of the 
dependent variable (u) or its derivatives. Since they are linear, any linearly independent com- 
bination of solutions will satisfy these equations due to the distributive property of linear oper- 
ators, i.e., 

L(xj +■ Xg + . . .) = L(xj_) + Uxg) + ... 


where 


L is a generalized linear operator 

This property may be applied to equation (17) as follows: Take Uj and u 2 to be independent 
solutions to equation (17). Next, define 

u 3 = a l u l + a 2 u 2 

where ai ui + u 2 = 0 if, and only if, aj = a 2 = 0, i.e., uj_ and U£ are linearly independent. 

Now substitute u^ into equation (17), then a 2 U 2 into equation (17), and add the two 
expressions [using (a) the shorthand operator 


dx* 


d2_ 

dy' 


d2 

. 2 
dz 


and (b) corresponding to the a^u^ solution and G 2 corresponding to the a2 U2 solution]: 


2 2 

V a i Uj + y a 2 u 2 + ^1 + ^2 ~ ^ 
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2 . 


Since v is a linear operator, and defining G = Gj + G 2 , we can write: 
0 

v (a^ Uj^ + ag u 2 ) + G ~ 0 


but 


a l U 1 + a 2 u 2 = u 3 


therefore. 


V Ug + G = 0 

Thus, u„ is also a solution to equation (17). Applying the same procedure to, say, equation (18b) 


d Uj 

d n 

b n 


K U 1 


h 

K U 2 


Adding 


Su l Su 2 h . , 

5 n + bn " K (u l 2* 


Since 3/dn is a linear operator, this can be written 

TiT (U 1 + u 2> - T (U 1 + u 2> 

But u 1 + u 2 = u 3 , then 
h 


Su 3 


K U 3 


Thus, U 3 also satisfies the boundary condition. 

B. ADDITIVE SOLUTIONS 

A useful ramification of the superposition principle lies in the fact that the solution to a 
complicated system may be formed by linear combinations of known solutions. 

Consider a rectangular plate with internal heat generation. The governing differential 
equation and boundary conditions are illustrated in the following sketch: 
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Explicitly, we have the system 
V 2 u(x, y) + G(x, y) = 0 
u(x,0) = f Q (x) 

u(x,a) = f a (x) < 19 > 

u (°,y) = g 0 (y) 

V u(l, y) = (h/K) u(l, y) 

The solution may now be written 

u(x, y) = UjU, y) + u 2 <x, y) + u^x, y) + u 4 (x, y) ( 2 °) 

The number of ancillary problems is taken equal to the number of nonhomogeneities in the sys- 
tem. Since the governing equation and the first three boundary conditions are not homogeneous, 
the number of ancillary problems is four. 

Substituting equation (20) into equation (19), we obtain the complete system: 


2 2 2 2 
V u i + V u 2 + v U 3 + V u 4 + g= 0 

u l +u 2 + u 3 +u 4 = f o 


( 21 ) 


u l +u 2 + u 3 +u 4 = t, 
“l + u 2 + “3 + U 4 = g o 


y Uj + v u 2 + 7 u 3 + V u 4 = ("j^) |^ u l + u 2 * u 3 + u 4 at x = 1, y - y 

The set of ancillary problems corresponding to this system can now be written as: 


Problem 1 


u^ + G = 0 


“1 = 0 
“2 = 0 
“3 = 0 


Problem 2 


v “2= 0 

u^ = 0 

u 2 (x, a) = f a (x) 
u 3 = 0 


Problem 3 


7 “3 " 0 
= 0 

u 2 = 0 


^(0, y)=g 0 (y) 


Problem 4 
V 2 u 4 = 0 
u x (x, 0)=f Q (x) 

U2 = 0 
u 3 = 0 


Vu 1 (l,y) = ^ u^y) vu 2 (l,y) =-^ u 2 (l,y) vu^l.y) = £ u 3 (l,y) yu 4 (l,y)=| u 4 (l,y) 


Note that each of the ancillary problems now contains only one nonhomogeneity, a much 
simpler form. Note also that their sum is equal to equations (21). 

Hopefully, we can solve each of the ancillary p 'Me ms or find the solutions in the liter- 
ature. Once we have these, we merely add them to obu*in the total solution to equations (19). 


It was mentioned that once a system has been degenerated to a set of ancillary problems, 
the final solution is the sum of the individual solutions. This holds, provided the ancillary 
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problems are properly defined. Care must be exercised in specifying boundary conditions so that 
the sum of the individual solutions equals the total solution. 

Consider, for example, the following three problems (taken from reference 10): 





CASE II 



Heat, which is generated uniformly over a circular disk S, spreads by conduction through 
a cylinder of height H and diameter D to a constant temperature heat sink. 

At first glance, the analyst might be inclined to assume that case III is the sum of cases I 
and II; he would be wrong. To understand why, we must correctly define the boundary conditions 
on each ancillary problem. 


The governing differential equation in each case will be the same 


S 2 u 


0 x' 


n2 

o u 

v 2 

by 


= 0 


( 22 ) 


The boundary conditions will be written in four parts, corresponding to the regions bl, b2, 
b3 and S . 


Case I: 


Uj(bl) = 
dUj (b2) 

^n 

duj (b3) 
a n 

&Uj (S) 

c>n 


0 

= 0 (adiabatic surface) 

= 0 

= Gj (constant flux) 


(23a) 


where n is an outward directed unit vector normal to the surface. 




Case II: 
du n (bl) 


0 


n 

Ujj(b2) = 0 

Sn (b3) 

“Tn = 0 

au II (s) 

Sn " G II 
Case in : 
u in (bl) = 0 

u in (b2) = 0 - • 

au TTT (b3) 

fn ■ 0 

Su m (S) 

an ~ °ra 

Now, add the governing differential equations for cases I and II 


aV 


a x* 


A 

a7 


a 2 u 


II 


ax" 


a 2 u 
ay 


ii 


= o 


Regrouping 


a 2 a 2 

_ ( U J + u n ) + (Uj + u n ) = o 

By our initial assumption Ujjj = u^ + u^. Thus, 


2 

a u 


m 


a 2 u 


III 


= 0 


a x " a y " 

The differential equation is satisfied. 

Next, add boundary conditions, beginning in region bl: 


Uj(bl) 


au n (bi) 


= o 


(23b) 


(23c) 


Note that we have mixed conditions which are inconsistent. This situation also exists in regions 
b2 and S (except that in region S we could define Gj + Gji AGqi to make the boundary conditions 
additive. Thus, case in is not the sum of cases I and II; we have, in fact, three nonanalogous 
systems . 
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C. FURTHER EXAMPLES 


An excellent set of examples of the application of superposition principles in the calculation 
of thermal spreading resistances may be found in Reference 1 in which the auther utilizes 
Green's function in the calculation of thermal spreading resistances. For a discussion of Green's 
function see Reference 11 . 
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SECTION VI 


SUMMARY OF OTHER TECHNIQUES FOR CALCULATING AND ESTIMATING THERMAL 

SPREADING RESISTANCE 


The two handiest ’’rules of thumb" relationships that can be used to calculate or estimate 
thermal spreading resistances were developed by Holm (Ref. 12) and Raillard (Ref. 13). Holm 
derives the equation for the thermal resistance of a circular isothermal source on the face of a 
semi -infinite slab as: 


R 


_J 

4 a k 


where 


(24) 


a = radius of circular source 
k = thermal conductivity of medium 

Figure 3 (Section II) shows the temperature profiles described by Holm's equation. Figure 
3 shows that 80 percent of the total resistance in a semi -infinite slab occurs within three radii 
of the source. It can be seen that, when the size of the source is small compared to the thickness 
of a slab of finite extent, Holm's equation can be used to make a conservative (high) estimate of 
thermal resistance. Such geometries occur often in microelectronic components. 

Raillard presents similar exact solutions for circular and rectangular sources having uni- 
form generation located on the faces of semi-infinite slabs. The equation for the thermal spread- 
ing resistance of the circular source bears a close resemblance to Holm’s equation: 


R = 


1 

7T a k 


(25) 


In Appendix B of his report, Raillard presents an exact closed-form solution for the rectan- 
gular source of uniform generation on a semi -infinite slab. He also derives the solution for uni- 
form generation in a strip of infinite length and finite width on a semi -infinite slab. 

Two references treat the case of the rectangular source of finite width and infinite length on 
a slab of finite depth. Wilcox (Ref. 14) treats the uniform heat generation source, while Gale 
(Ref. 15) presents thermal spreading resistances for the uniform temperature source. The re- 
sults of these two studies can also be found in Reference 16. 

Muller (Ref. 17) and Kennedy (Ref. 10) present exact solutions to the problem of a circular 
source at one end of a right cylinder with conduction to the side of the cylinder, the other end of 
the cylinder, or to both places. Kennedy’s source is uniform while the intensity of Muller’s 
circular heat source varies exponentially within the source region. 

Hein (Ref. 1) examines stady -state heat transfer in a rectangular substrate or slab having 
multiple heat sources. He has integrated circuits in mind. He considers convective heat transfer 
and lead conduction for various heat-sinking conditions and his solutions are mathematically 
exact. 
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Finally, it might be well to point out a principle which is somewhat analogous to Saint - 
Venant’s Principle in the theory of elasticity. Simply stated, T, the temperature in any thermal 
spreading resistance problem, varies with 1/L where L is the distance from the source, when 
L is large compared with the characteristic demension of the source. That is, 


T oc 


J. 

L 


when L » a where a is characteristic source dimension 


This is indicated by the form of the solution of Fourier's equation for spherical flow, i.e. , 


q 


A dt 
dt 


«i - 1.1 


f 


2 ^ 
r dr 

. 2 
Attt 


4;rk m (t l -v 

04--# 
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SECTION vn 


PLANNED APPLICATION OF COMPUTATIONAL TECHNIQUE 


General Electric’s Aerospace Electronic Systems Department has developed an extremely 
compact and thermally efficient packaging configuration for computer circuitry. Integrated cir- 
cuits are bonded to multilayered printed wiring boards which are in turn bonded to compact 
forced-air-cooled heat exchangers. This configuration is illustrated in Figures 18 and 19. 

A heat transfer analysis program that calculates the temperature of each flatpack using 
Gauss- Seidel iteration is currently employed. This technique has proven costly: 200 + iterations 
are required. Computer costs of $25 to $50 per analysis have been experienced. 

The exact computational technique described in this report will be implemented in the near 
future for this type of analysis. Costs are expected to be an order of magnitude lower than those 
with the iterative technique. (See Figure 9 of Section III. ) 
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Nl = 16 



Figure 18. Memory Board Module 


2ND ROW 
(1 = 3) 


1ST ROW 
0 = 2 ) 


INTEGRATED CIRCUITS 



Figure 19. Location of Conductances and Temperatures 
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SECTION vm 


RECOMMENDATIONS FOR FUTURE WORK 


1. The program should be rewritten for three-dimensional fields. Special attention to 

such techniques as using the method developed in this report for the main coefficient 
matrix on the submatrices themselves should be examined. The fact that very large 
matrices can be efficiently handled by this technique should be capitalized upon. 

2. The program as it now stands should be used in a thermal spreading resistance 
parametric study similar to but larger in scope than the illustrative problem of this 
study. Nine hundred to sixteen hundred finite elements should be used in such data 
generation. 
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APPENDIX 


GENERAL FORTRAN Y VERSION OF COMPUTER PROGRAM 
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% 

$ 

$ 

JDENT 727-9C4* KELLY NEC ,D6L-B ' 

OPTION FORTRAN 
—USE MEMORY/1QOG/ 

65078300044700 

I 

% 

I 

entry main 
FORTY 

IWCOnF I0MF 


CHAIN 

MAIN 

SUBROUTINE MAIN ' 
_PARAKFTPR MAXCORslflnn 



COMMON /MEMORY/CORE (MAXCOR) 



COMMON /FILES/ INFILE, 10FILE, I FI LEI# IF I LE2 
-COMMON /TlMES/lTlME(91.inATFf2).f!EI TtM 


COHMCN /DEBUG/ I DEBUG 
PARAMETER MAXMAT*69 

&A-B.AMF TF R MAXHFrs* ... 


DIMENSION lOFF(MAXOFF) 

DIMENSION KARDU4) 

£A&A££XEB_ H A X T I T * 1 2 

DIMENSION I TITLE ( MAXT I T ) 
PARAMETER MAXTIM«2 

BATA KPEBUQ/5HQF8Ufi/ 

data iqebug/i/ 

DATA INFILE/05/ 

... DATA I QFILE/ Q hi 

DATA I F I LEI/ 0 7/ 

DATA JFILE1/6HU00007/ 

DATA IF1LFP/DB/ 

DATA ISIZE/4HSIZE/ 

DATA ICORE/HAXUOR/ 

DA TA KTIT L E/6HTIT LE / 

CALL FXOPT< 67 # i# 1< t) ) 

CALL FXOPT(68,0,U # Q ) 

C-AL.L-.F XQP T (Ag-iJ-iJ-i-Q..) 

CALC FXQPT ( 7 0 » d » 0 # 0 ) 

CALL FXOPT ( 71 # U # D # 0 ) 


IF(KARDCl).NE.KTITLE) GO TO 94 
CALL SUPERT<HAXTIT,KARD<2>> 

SLl-gM.Il?iUE 

call newlin 

SPRITE! IDT ILEi 7d)KARD 
CALL TIHDATCITIME, I DA TE > 

“ CALL EL T I ME ("OELfTH >’~ 
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.c 

c 

c 

JL 

c 

c 

JE. 


5S8 continue 

REAO( 1NFILE # 77JKARD 

_XLX!JR M A T<1 3 A 6 , A 2 ) 

CALI NEWLIN 
WRITEU0FILE*7a)KARD 

— 78 _F0ftriAT'( 2H ♦ , 13A6, A2, 1H* ) 

IF(KARD(1KNE*KDEBU0) 00 TO 557 
IDEBU0SO 
CO TO 558 

557 CONTINUE I 

DECODE(KARD,l)KEr,N 

1. FORMA T( A6, 4X, 15) _ 

IF (KEY. Eii. ISIZE) 00 TO 2 
CALL NEWLIN 
WRITE ! 10F1LE,55) 

55 FORMAT(40H THE ABOVE CARD SHOULD BE A ’SIZE 1 CARD.) 

STOP 

2_ IF(N*LTtMAXHAT) 00 TO 222 

CAUTNEWLIN 

WRITE ( I OF I LE* 223 ) MAXMAT 

223 FORM A T ( 18h SIZE QREATER THAN, 110) 

STOP 

222 IF(N.GT.O) GO TO 224 

CA LLJLEHL I N 

HR1TE{ I OF I LE, 225 ) 

225 FORMAT (17H SIZE LESS THAN 1) 

STOP 

224 NSQsMN 

CALL SETNUM!2*N*3) 

RECO RD S ARE IN SYST EM STANDARD RANDOM FORMAT 

WHICH MEANS THAT IF A RECORD IS GREATER THAN 318 WORDS 

THEN THE RECORD WILL BEGIN IN A NEW BLOCK 
UUL-EMP JL.alO.CK EVERY T I ME 

I TRIED USING PURE DATA RANDOM FILESC11-13.72 ) , BUT 

FOR, SOME REASON They did NOT SEEM TO WORK PROPERLY. 

*-8L0CK=< (NSQ-D/318 ♦ 1)*(2.N ♦ 3) 

CALL NEWLIN 

*riteuofile,8;)nblock 

87 FORMAT! 110, 43H BLOCKS OF RANDOM DISC STORAGE ARE REQUIRED) 
tiUNKS*<NBL0CK-l>/12 ♦ 2 
CALL NEWLIN 

W R I TElTOF I LE, 8 6 )NLINKS “ " 
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86 

ai 

FORHATC I10,42H LINKS OF RANDOM OISC STORAGE ARE REQUIRED) 
NTlHESaO 

-CALL ftp THQR ( 1 , I ERR, NL 1 NKS. JF I LEI ) _ __ 



c 

If ( IERR.EQ.0 ) GO TO 88 
request was REFUSED 
-CALL__tiEWLlN . .. 



82 

«R I TE C I OF ILE# 82 ) 

F QRMA T ( 2 9 H REQUEST FOR DISC HAS REFUSED ) 

JiTT^P<gWTtHCC*1 



IF<NT!«ES.LT.HAXTIH> GO TO 81 
CALL NEWLIN 

StBJ-tE! I0FIL£,83) NT1HES 

83 

F0RMATC25H REQUEST FOR DISC REFUSED# J8#16H TIMES. GIVE UP.) 
STOP 



dfi- 

-C-Ai L KFTSTZfMSQ) 



c 

SEE IF WE HAVE ENOUGH CORE 
MC0RE=5*NSQ ♦ n 

CALL NEWLIN f 



1Q1 

WRITE! IOF1LE, 101 ) ICORE 

F0RMATU1D,38H WORDS OF CORE ARE CURRENTLY AVAILABLE) 
CALL NEWi. IN 



102 

WRITE! I OF ILE# 102 JMCORE 

FORMAT! 110, 39H WORDS QF CORE ARE REQUIRED FOR THE J08 > 
Jf.IM.QR E jJ. EPICURE) GO TO 3 



C 

c 

c 

GET MORE CORE 

ICORE IS AUTOMATICALLY UPDATED TO REFLECT THE ACTUAL NUMBER 
QF WORDS THERE ARE UPON RETURN 



c 

c 

THE ROUTINE DOES NOT FAIL 

IT KEEPS TRYING UNTIL IT GETS THE CORE IT HANTS 
CALL NEWLIN 



103 

WRITE! IOFILE, 103) 

FGRMAT(33H GET THE ADDITIONAL CORE REQUIRED) 
CALL GIMMF(MCQKE, ICORE. CORE} 



C 

3 

COMPUTE LINEAR OFFSET FOR EACH MATRIX 
I OFF ( 1 ) =1 

-0.0-. gi-MX of f„ 



4 

C 

I OFF! 1 >* I OFF C !-l)*NSQ 

Ni)N*N 

RESERVE SPACE FOR 2*N MATRICES 



91 

JsaVNUM r ~~ — 

DO 91 1 = 1, J 
call NEWNUM(K) . 



call FIRST(NUM#CORE(ICFF!1)),NUH,NUM,CORE(IOFF(2) ) # NUM# HUM, 
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CALL FINISH 


CORE (I OFF <3 >>, HUM# NUM# CORE < |0FF(4 )) * NUN, NUM# 
C0ftEUCFF(5)),NUM,NU«#C0RE(I0FF(6))#NUM) 


STOP 

end 

-CPiRSJ__ F1RST„ ; 

SUBROUTINE FIRST! NUM# 01/ NRl, NCI, 02, WR2, NC2, 03, «R3, NC3/Q4, NR4, NC4# 
1 Q5,NR5 ,nC5, 1VEC/NVEC) 



“Common /DTb u gTTdFbUS 

DIMENSION IVEC(NVEC) 
REAL 01(NR1,NC1) 



REAL Q2(NR2,,s ; C2) 
REAL 03(NR3,NCJ) 



REAL Q4(NR4 # NC4) 



REAL 05<NR5,KC5) 
DIMENSION I Cl (4) 



DIMENSION 102 (4) 



DIMENSION 103(4) 
DIMENSION 104(4) 



DIMENSION 105(4) 



iai(l)aO 
I G2 ( 1 ) = 0 
J Q.3iU__=_Q_ 


I 0 4 ( 1 ) s 0 
I 05 ( 1 > = 0 
-■IC1(2?«0 


102(2)30 

103(2>o0 

lii<2J?0 
1 05 (2 )«0 
I01(3)«NR1 
I Q 2 ( 3 ) * N R 2 


IJ3(3)sNR3 
IQ4(3)=NR4 
115X3 > = NR_5 

101(4 )«NC1 
102(4 >«wc2 
103(4 )SNC3 


104(4 ) = NC4 
105(4 ) = NC5 

c A lk S£ TUPA ( 0 1/ N R 1 , N C 1 / C2 / r»R2,NC2 . 0 3 / *R3,M C3, } 4 ,_Ng4j tl C 4 /_ 
1 Q5,NR5,nC5 # IV6C/NVEC/NUK) 


( 
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CALL HSUB<Q5, I-J5# 04 , IQ4#Q4# IQ4* IERR) 
CALL MINV{ Q4, NUM,NUM# I VEC, OET) 


CALL W«PY(Q4 # ]U4#05 # IQ5,Q2» IQ2, 1ERR> 
CaLL SAVE<Q2, 1 02, NAHE + l) 

r A t i kMDY/n-* t n -» n- , ~ - 







CALL RSTP< <J2, 132, name-1 ) 

__«2 ■_.?_# . 

call MADD<Q1, IQliOZ# 102,02, IQ2/IERR) 
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ofo o 


IFUERR.NE.OJ CALL 0 U I T C IERR) 
CALL SAVE(Q2, 132, NAME-1) 

81 ; 

IF<NAME,GE.4 ) GO TO 20 
IF< ICEBUG.NE.O) 90 TO 667 

C ALL NEWLI N 

WRITE ( 10FIL6.567) 

567 FORMA T (14H DEBUG PQfST ’89 

CALL HATHAT 

667 CONTINUE 


PRINT OUt THE RESULTS " 

CALL NE WPAG 

CALL NEWLIN 

NRITE(IOFILE,8801) 

6801 FORMA T (9H A MATRIX) 

CALL R STR (Q l, IQ1,1) 

CALL PMATCOl, I Ol (1 ) , 101(2)# I Ql ( 3 )', 131(4)} 
CALL NEWPAG 



CSETUPA 


RETURN 

end first 

ENn 


SETUPA 


SUBROUTINE SETUPAOUNRH# NCH# V, NRV# NCV# C# NRC, NCC# X, NRX, NCX, 

■i— - - Y«RRY iNCY i I VEC f HVECj,NUH? 

COMMON /FU.ES/lNFILE#IOFILE# IFILEl# IFILE2 


COMMON /DEBUG/ UE8UG 
■HiMENSION _l VEc YNVEC ) 
DIMENSION H< NRH# NCH) 
DIMENSION V(NRV#NCV> 
JUM£ Mli)N--C{NRC,NCC> 
DIMENSION X<NRX,NCX) 


dimension Y<NRY,NCY) 


■aiMFNSION I M C 4 ) 

DIMENSION I V ( 4 } 
DIMENSION I C ( 4 ) 

DIMENSION I X ( 4 ) 

DIMENSION ! Y ( 4 ) 
DIMENSION KARD114 ) 
SATA 1HEAT/4HHEAT/ 
IH(1 )»NUM 
IH(2)=NUM 


1H <3 )sNRH 
I H ( 4 > = NCH 


tV<l)=NUh 
I V ( 2 ) a NUM 
IV(3)sNRV 
I V { 4 ) *NC Y 
icq )»NUM 
I C < 2 ) a 1 


IC(3>«NRC 

IC(*)=NCC 

IX(1)=NUM 
I X ( 2 ) = NUM 

I X ( 3 ) sN RX 

IX<4 ) = NCX 
I Y< 1 )*NUM 



I V < 3 )*NRY 
I Y ( 4 ) aNC Y 

CALL NEHNUM(1M1/£X.CJ_... 
CALL NENNUM< INOEXH) 
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CALL NEWNUH(INDEXV) 
CALL MATZER<H,IH> 
CALL MATZER(V/IV) 


^ALL MATZER < C , IC > 
PI-3. 14159265 
C 1 = • 5 / P I 


R1sC1*(AL0G(2.0>)/(2,0#P!) 
R2bR1/2* 0 * “ 


H<NUH,1)*1,Q/R1 


IBP 

DO 180 H»2# MUH*1 
H(jui>32,o«Haa) 


210 

DO 210 M*l# NUM 
V(H,l)aPI/2 f 0 
DO 310 Hal i NUM 


DO 300 N*2,NU« 

H<M*N)«l,fl/Cl*e/(4,9*M>**LOO<fLOAT(N)/<rLOAT(H>^l*0) )) 
YLBr N>»P 1*4. D«< FLOAT 0 ) " * 

23JL 

IF<H.OT,l) GO TO 2BI 

H(l,K)5 t 5#H(l,»0 

h<Nb^.N)».5#HlNUM»N) 


30 0 

m_ 

V(M,lO>=Pl*<4.Q* FLOAT ( N >-5. 0 I/Z.O 

continue 

CONTINUE 



DO 340 M=1,NUM 


340 

H ( H# NUH)=0«0 
DO 370 n»i;num 


370 

v ( 1 , N ) » Q . o 
C2 a l • 0 
C3sl , 0 


DO 460 M=1,NUM 
DO 450 N=1,NUH 
V(H#N)=C3*V(M,N> 

450 

460 

H(H#N)=C2*H(M # N) 

CONTINUE 

CONTINUE 


461 

READ ( I NF | LE, 461# £nD« 468 >KARfl 
FORMA T( 13 A6 # A2 ) 

Call newlin _ 


464 

HRITE<I0FILE*464)KARD 
FORh'AT(2H * j 13A6# A2 * 1H* ) 

. I FCK A R D 1 ) *E O.I . HE AT). GO TO 462 - 


CALL NEUL1N 





59 



WRITEU0fILE,463> 

463 FORMAT (35H ABOVE CARD SHOULD BE A * HEAT* CARD) 
-§-T0p 


462 


D6C0DE(KARD, 465)11, J2 
465 F ORMA T (10X#2I5) 

lFjJJL t_LE t QtOR. U,ST, !2 ) QO TO 479 

I to nr tiriut t>n *n jt-i 


IF < I2.GT.NUM) 00 TO 473 
GO TO 479 
_472 CAl I NFVI 


I TE C IOFILE# 474 > 

474 FORMAT (20H II IS OUT OF BOUNDS) 

Stop 


473 CALL NEWLIN 

HR I TE( I OF I LE# 4 75 ) 

475 FfiRH A TfPnw t? 14 nut nr ammpg) 
STOP 

479 CONTINUE 

— Q-0— 4 6 6 - 

466 C(N,X)*0*0 

DO 467 N» II# 12 

-167 C(N#1 i*V(2#N)»2>n - 


CALL NEWPAG 
URiTFMnrn N c r 4ft<n 


469 F0RMAT<26H INITIAL HEAT INPUT MATRIX) 
CALL PMAT<C#lC(l>#IC<2)#ICC3),IC<4)> 

-Q.Q Tn 4? ) 


468 CALL NEWLIN 

470 FOR*ATU9H UNEXPECTED END OF FILE* 
WRI LE < I OF I L E« 4 7 0 )_ 


EXPECTING A »HEAT« CARD) 


STOP 
471 CONTINUE 

CALL NEHPAQ 


CALL NEWLIN 
MR I TE C I OF 1 LE, 701 > 
_70.1 F0RMAT(9H H MATRIX! 


CALL PMAT(H,1H(1),1H(2), I H ( 3 ) # I H { 4 ) ) 
CAU^NEWPAG 


WRITE(I0FILE,702) 
702 FORMAT! 9H V MATRIX) 

call 


CALL S A V E ( C # I C # I N D £ X C ) 
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CALL SAVE(H # IH, INDEXH) : 

CALL SAVECV, IV* INDEXV) 

_J/<IOEBUO.NF.D) 60 TO 664 

CALL NE ML IN : 

NRlTE(lOF!LEf564> 

-£0 RHAT(14H DEBUS POINT A) 

CALL MATKAT 

CONTINUE 

"The >ATR ,C ES ON THE DIA0ONAL (AND OF COURSE THEt« INVERSES) 
So « B f N L,2 NLY F0R THE rlRST ' SECOND, AND CAST TIHES. 

IQ0«1 “ — 

IF(H,E0.1.0R.N.E8.2.0R»M,60,NUN) I80«e 
IX<n»NUH 

IX ( 2 ) * N U M — 

CALL MAT2ER(X # IX) 

I Y ( 1 ) bNUM 

I Y ( 2 ) a ,Sj U H 

CALL HA?ZER<Y,JY> 
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CALL MMPY<C,IC#X#IX,y, IH# IERR) 
IF(IERR.NE.O) CALL QUIT(IERR) 
CALL HNEGIH, I H > 


C 

SAVE H (B#D#F##.«) 
call SA VE (H# I H# M*2 ) 
—If (M.fiT.l) ftn T 0 .... 888 


CALL RSTR ( H# IH# I NDEXC ) 

CALL MMPYIC# ICs H# IH»X# IX, IERR) 
If < IERR.QT.O) CALL OUlTMERR) ... 


C 

SAVE A 

CALL SAVE ( X# I X# M*2*l) 
GO TO 1150 



888 

CALL MMPYIC# IC#Y#IY#X# IX# IERR) 

call mnegix# I X ) 

IF (I ERR. NE,_0 3-CALL QUIT! IERR) 



HI 

SAVE (C#E#G# ... ) 

CALL SAVE<X, IX#M*2-1) 
CONTINUE 


IF ( I DEBUG* NE • 0 ) GO TO 665 
CALL NEHLIN 

WRITE(I0FILE,S65) .... 


565 

FORMAT (14H DEBUG POINT 8) 



$A5 

CALL MATMAT 

CONTINUE 



C 

RETURN 
END SETUPA 

_£M 



csave SAye 

SUBROUTINE SAVE< A# IA# INDEX > 

COMMON /FILES/ INFILE # I OF HE, IFILE1# 1FILE2 

DIMENSION Ad) 

DIMENSION I A ( 4 I 

pjii.gj.!sj.. i o,N lam ; 


DATA NRW/0/ 

DATA kount/o/ 

I F( IND EX . LT.l. Q R.IND E X . GT. XQUNT) GO TO 801 


C 

IFC IA(1),LE.C ) CALL QUITd) 

IF( IA(2).L6.0 ) CALL QUITd) 

CHECK TO SEE IF RANDOM RECORD SIZE IS LARGE ENOUGH 



IF((IAC1)«IA(2)).GT.MAXWRD) CALL QUITd) 
I B d ) = I A d ) 

I B ( 2 ) a I A ( 2 ) 



I B {3 ) = 1 
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HRITEUFILE2UNDEX) IB 
NR-IAU) 

NC»IA(2) 

RCOLSMAtA) 

I3«NCQLS*(NR«»1)*1 _ 

WRITE? IFILE1 » JNDEXXUt !)# J*U# II* 12 ># Ilsl# 13 # NCOlS) 

NftWaNRW*l. 

RETURN 

c 

ENTRY RSTRU#IA, INDEX) 

IF < IND E X.LT.l .QR. IND EX.G T » KOUNT ) 80 TO 601 

READ? IF ILE2 U NDExTIB 
IF(IB(3J*EQ«0) CALL QUIT(INBEX) 

C CHECK TO SEE IF THE DIMENSIONS OF THE MATRIX ARE LARGE ENOUBH 

IF(I0(1) # 0T*IA(3>) CALL OUI TCI > 

IF ( IB(2).GT* I A(4) ) CALL OUITU) 

I A ( 1 ) a IB < 1 ) 

I A { 2 ) * 1 8 ( 2 ) 
f-’R« I A Cl ) 

hSLRlMZ) : 

•vCOLSa I A ( 4 ) 

I2*NR-1 

13 = (<CQLS*( N R - 1 ) » 1 

READ? IFILEIMNBEXX ( A < !> # 1 « U# II* 1 2 ) # ll»l# 1 3 ♦ NCOLS > 

NRW=NRW*1 

RETURN ' 

C 

ENTRY FINISH 

CALL NEWPA O 

CALL NEWLIN 
NR1TE(I0FILE#U11)NRW 

nil FORMAT ( 1 1 0 > 21H READ-WRITES EXECUTED? 

CALL NEHLIN 
WRITE(1QFILE,6> 

6 F0RNATC12H NORMAL HALT) 

STOP 

C 

ENTRY SET 5 I Z ( NUM8R ) 

HAXHRO=NUMBR 

CALL RANSIZ(IFILEl#MAXWRD#0> 

CALL RANSIZ(IFILE2,3) : 

RETURN 



c 


ENTRY SETNUM(NUMBR) 
MAXKNTsNUMBR 


c_ 


CALL SETDUM(NUNBR) 

return 


ENTRY NEkNUMf NUHBR) 

kount*kount*i ■ 

—IF t KPllMT.LE. MAXKNT1 Eft TO 5f)l 


502 

wRlTEUOf ILE,502)MAXKNT 
FORMATUOH MORE THAN, 110, MATRICES) 

-CALL QUIT ( 1 ) 



501 

CONTINUE 
NUMBRaKOUNT 
-1B<1) = 0 


1B(2)=0 

IB<3)=0 

wRITEtinLE2»KOUNT)IB 



RETURN 







1- 

NSOa I A (3 ) • I A ( 4 ) 
00 1 in, NSO 
iU.) = 0,0 




return 


C_ 




ENTRY NNEG ( A # I A > 
KSQ* I A (3 )• l A ( 4 ) 
DO 82 1=1, NSO 

C 

82 

AU)=-A<I) 

RETURN 


ENTRY MAT IDN( A, I A, NUM ) 

IFUA(3).LT,NUM.0R,lA(4),LT t NUM) CALL QUIT<1> 
IAL1) = NUH 

I A < 2 ) = NUM 
K'SQ= I A ( 3 ) • I A ( 4 ) 

t?Q 2 1=1, NSQ 


2 

A ( I ) =0 • 0 
DO 3 1=1, NUM 
J'«SQ=IA(3J.(1.1) ♦ I 



3 

A (NSO ) =1 , 0 
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RETURN 

c 

_ 6Q1 C ALL NEHLIN 

" RITE(10F[LE*802)INDEX 
602 FORHAT(ZOH BAD MATRIX INDEX OF#UO> 

CALI QUIT(1) 

STOP 

C END SAVE 
END 

CNEWUN nehlin " 

SUBROUTINE NEHLIN 

COMM ON /DATTIM/ITIHE(2), I DATE <2 ) 

PARAMETER MAX»20 
DIMENSION lHEAD(MAX) 

DIMENSION JHEAD(MaX) 

DIMENSION IVECUiHORDS) 

DATA JHEAD/HAX*6H / 

DA TA IOFILE/6/ ~ 

DATA MAXLIN/55/ 

DATA I NAME/1H / 

— DAlIA NL 1 NES/fl/ - 

DATA IG01/I/ 

DATA IPAGE/0/ 

am.jLZ/t* 

DATA IBLANK/1H / 

1002=1 

0 0-10 < l£32?/IQ01 

32 IF<NHNES.GT.O ) QO TO 5 
I PAGE* IPA0E*1 

WRITEdOF n, E * IAJ 

18 FORMAT (1H1# // ) " " : 

HRITE<I0FILE,6) IT1ME# IDATE.IPA6E 
fl FORMAT! 9 H TIME IS ,2A6 , 4H ON , 2A6,80X. 

1 5HPAG£=, 15) 

WRITE! I0FILE,87)JHEAD 

87 FORMAT! 1 0 H . TITLE*, 2Q A6) 

HR I TE < I OF I LE# 7 ) I HEAD 

7 FORMAT! ID H SUBTITLE*, 2QA6) 

DO 8 I«1,N2 

8 WRITE! IOFILE# 9) 

9 FORMATUH ) 

NL.1 NES»N 2j l 3 

5 NLINES=NLINESd 
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s ' 

IF< N LINES.G7,MAXLIN)NLINESifl 

return 


entry newpag 

NL!NES=0 

RETURN 

c 

entry sethedcmhoros* IVEC) 

-1Q02b? 


1 

go TO < 1 , 2 > # 1001 
CALL TIMDAT( l T I ME, I DATE > 
_I GDI s2 


2 

ML1NES=Q 

IPAGEaO 

DO 3 Id .MAX 



IwKllinf 


4 

fi 

IHEADC I >«IVEC< 1 ) 
RETURN • 


80 

ENTRY SUPERHNWORDS, IVEC) 

DO 80 Icl,NWOR0S 

-JH£ADm«lVECIl> 


c 

RETURN 
END NEWLIN 

JLlifl 


CQUIT 

QUIT 

SUBROUTINE QUITUERR) 

COMMON /FILES/ INFILE, I OF 1 LE, IFILEl, IFILE2 


DIMENSION 18(3} 
NSOao 

NSObnSQ**NSQ 

21 

WRITE! I OF I LE, 21 > I ERR 

F0RMAT(22H QUIT BECAUSE OF IERRa,UO> 



GO TO 777 


c 



777 

1 60 = 0 
CONTINUE 

-Ng.lIEJ LQf.l.LfLil? ) 


12 

FORMAT(19H MATRIX INFORMATION) 
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WR1TE(I0FILE,13> 

13 FORMATOOH MATRIX NUMBER NO. ROWS IN US E NO. 

1 USE EVER USED) 

COLS IN 

uu 10 * S 1 , MA XKN T : 

READ ( IFIL62 1 1 JIB 

. i° wf?iTeuorae,n>i , ib 

11 FORMA TM 1 20 ) ' 

IF ( I G0» EQ » 0 ) RETURN 
CALL PDUMP 

c 

& I OP 



ENTRY SETDiWUERR ) 


C 

HAXKNT*IERR 

return 
End oy 1 7 


fcND — 

•INVRS SUBROUTINE TO OBTAIN INVERSE OF MATRIX, CALL OECOM 

- CD600D4.Q07 DATE 05/04 /«* 

FIRST INVRS020 
INVRS030 

• 

iUBKUUlINfc !NVKS(A,INTR,MSIZE,NNJ 
DIMENSION AIMSiZE.MSIZE), INTRIMS I ZE> 

Srjgfciax ^ V f RSE 0F PECOHPOSED MATRIX 

IVVRS040 
I NVRS05 0 
IMVRS060 


subroutine decom must be Called First 

NcNN 

IF(INTR(N))18»17,18 

INVRS070 

INVRSOBO 

InVRSobo 


1 DO 13 K*1,N 
KM«K-1 
IF < K M ) 2. 7. 2 

INVRS1D0 

INVRS110 

F N VR<Si •> n 

• 

2 IF(K-N)3,7,3 

COMPLETE reduction below diagonal 

3 KP»K*1 

INVRS130 
INVRS140 
I NVRS150 


DO 6 1=KP,N 
X s A ( I , X ) 
IF(X)4,6,4 

INVRS160 

INVRS170 

INVRS180 


4 DO 5 J*l, KM 

5 A( I # J)3A(K»J)*X*A( I , J ) 

6 CONTINUE 

INVRSi«g 
INVRS200 
INVRS71 n 

• 

DIVIDE THROUGH BY PIVOT ELEMENT 
7 X*1.0/A(K,K> 

4<K,K)»1. D 

I NVRS22 0 
INVRS23D 


DO a Jsl,N 
8 A(K,J)=A<K, J)*X 
_ IF (K.M >9.13,9 

INVRS250 

INVRS260 

• 

REDUCE TERMS ABOVE DIAGONAL 

I N V R S 2 8 :1 
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9 

Do 12 1 s 1 , K H 
X*-A(I,K) 

— LElLXlltJ 10 




INVRS290 
INVRS3C0 
1 NVRS31 n 

10 

U_ 

A(I#K>80.0 
DO 11 j«l,N 

JLLlt i , j ) 



INVRS320 

INVRS330 

TNVRS34A 

12 

13 

• 

CONTINUE 

CONTINUE 

_IN TPRPM A N RE rni HUNS 



INVRS350 

INVRS360 

1 NVRStWn 

• 

KM=N-1 FROM PREVIOUS LOOP 



INVRS380 


DO 16 J=1,KM 
JLoW-.l 




I N VRS390 

TNVRSAdn 

U. 

KPbINTR(K) 

!F<XP)14,16,14 

nn i * i si f n 




INVRS410 

INVRS420 

l M V R S A .1 n 

I5_ 

X«A(I,KP) 

ACI,KP)sAU,K) 

.A(I,K)=X 




INVRS440 
INVRS45 0 

I N VKS460 

16 

17 

-16 

CONTINUE 

return 

KMeN-1 




INVRS470 

INVRS4S0 

INVRS490 


DO 21 I *2 , KM 

KsiNTR(I) 

IF<K)l9,21,19 




INVRS500 

1NVRS510 

INVRS520 

19 

KPM-1 

DO 20 Jsl , kp 
_X* A i I » J ) 




IN VRS530 
INVRS540 
I N V R S 5 5 fl 

20 

21 

A(I,J)*A(K,J) 
A<K, J>*X 

continue 




INVRS56Q 

INVRS570 

INVRS580 


80 TO 1 




INVRS590 

mVHI 

end 



ROUTINE 

1NVRS600 

MINV0fl2n 

« 

subroutine MINVU 

--DIMERS ION ACMSIZE 

CD600D4.007 
,MSIZE,N, INTR,DET> 

DATE 05/04/65 

MNVQ030 
fl 1 N V 0 0 4 0 
m jNVflnsn 


CALL DECOMU* INTR 
CALL DTMNU, INTR, 
-CALL INVRSU, INTR 

#MSIZE,N> 

MS I Z£, N, DET > 
^MS_IZE,N> 



M I N V 0 0 6 0 
MINV0070 
ti I N V 0 0 R 0 

•MKPY 

RETURN 

END 

CD600D4, 

(U2 MATRIX 

MULTIPLY 

ROUTINE 05/18/66 

M I N V0 130 
M|NV0140 
HMPY0C02 

• 

COPYRIGHT 1966 BY GENERAL ELECTRIC COMPANY 

MMPY0Q03 
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SUBROUTINE MMPY(A,iDA,»,ID9,C,IDC,IND> 

DIMENSION AU>,BU),C<1),IDaM),ID8<4>,IDC<4> 
!A«!DA(1> 

JAriDA(2) 

HA«IDA(3) 

jalTj db (2 ) ~ 

MB«IDB(3) 

LPs I DC ( 1 ) 

JC“ I D C ( 2 ) ~ 

NCb I DC (3 > 

*C|lI PC(4) ' 

I ND“0 

IF< JA, NE# 18 M Nflo2 

Lf UM.LT.U)»0P.(NC.LT.J8niND»l 

If ( 1 ND« NE • 0 ) RETURN 
IDC(1)=IDA(1) 

I PC < 2 ) = r DB ( 2 ) 

DO 1 1*1, IA 
DO 1 K*1 » J8 

; 

C(LC)*0,0 

K-l > 



LA=MA»(J“1)*I 

LB=KB+J 

1 C(LC)=C(LC)*A(LA>*8(LB) 

return 

END 

• HADP CD60004,010 MATRIX ADD ROUTINE 05/18/66 

• ~COPYRTQfiT '1966 0Y GENERAL ELECTRIC COMPANY 

SUBROUTINE MA DD ( A , I DA . 8 , l DB, C , 1 DC, I NO ) 

DIMENSION A(1?,B(1),C(1),I0A(4?,1DB(4),I0C(4) 

I A ■ ! T A ( 1 ) 

UA= | p A ( 2 > 

18 =108(1) 

JB=!D8(2) 

ic=ioc<i) 

JC* IOC ( 2 ) 

MA* 1 D A ( 3 ) 

MB a I [ 8 ( 3 ) 

MCc I DC (3 ) 

NCs J DC (4 ) * 7 


HMPY0040 
MMPY0050 
MMPYD060 
MMPYO 070 
MMPY0080 
MHPY0090 
MMPY0100 
MHPY0110 
MMPYO 120 
MMPYQ130 
MMPY0140 
MMPY0150 
MMPY0160 

MNP Y 018 0 
MMPY0190 


MMPY0203 

AHPY 022 0 
MMPY0230 
MMPY0240 


MMPYO 26 0 
MMPY0270 
MHPY0280 
MMPY0290 
MMPY0300 
MADD0062 
MAOOOOG3 
MADDO 0 40 
MADD0050 
MADDO 060 
MADD0070 
MADDO 080 
MADDO 090 
MADD0093 
MADDO 096 
MADD0100 
MADD0110 
MADD012 0 
MADD0130 
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r~ ■' 





I nd»o 

1F < ( IA.NE. 18). OR* ( JA.NE, JB) ) JNDs2 
— LE-LlHC# LT #IA).OR«fHC*LT»JA))IND»l 

MADD0140 

MADDMAiT 



IF(IND.N6,0>RETURN 

IDC{1)=!DA(1) 

— LHC ( 2 >= I D A ( 2 ) 

MADD0170 



DO 1 J«l# JA : ; “ 

-Lftpf lAtJH 

MADD0180 

MADD0190 



LBbMB# JM 
LC=HC*JM 
-&°-l I=iiIA 

MADD0210 

MADD0220 



LAsLA* 1 
LB»L8*1 
— Lc si r*i . 

HADD0240 
M ADD 025 0 


1 

C<LC)*A<LA>*8<LB> " " “ — 

return 

_6Nn 

MADDD270 
MADD0280 
MAi>nn ?on 

MSUB 

* 

CD60004 t «l MATRIX SUBTRACT ROUTINE 05/18/66 

COPYRIGHT 1966 BY GENERAL ELECTRIC COMPANY 
-JjU.3 ROUTINE MSUB f A* TDAj9« IOB.Cj I DC » INDi 

MSU80002 
MSUB 0 0 03 



DIMENSION A<1),B<1),CU),IDA(4),IDS(4),IDC(4) 
IAsIDA(I) 

_*L4 = IDA(2) 

MSU8O05O 
MSUB 0 06 0 
M s UH n ft 7 ft 



IBsIOBU) ~ 

JBs I DB ( 2 ) 

_LC« I DC < 1 ) 

MSU80060 

MSUB0090 



JC« 1 DC ( 2 ) 
MAs I DA ( 3 ) 
MBs 108 ( 3 ) 

MS U 80096 
MSUB01D0 
MS Lift ni i n 



MC ■ I DC ( 3 ) 
f ,i C«IDC<4) 
-LaO.flfl 

MSUB0120 

MSU0O13O 



IF( ( IA.NE. IB). OR* ( JA.NE. JB ) ) I NO -2 
IF< (MC.LT. I A ).0R. (NC*LT, JA) > IND*1 

MSU80160 

MSlIflm 70 

I DC < 1 > a 1 DA C 1 ) -* 

I0C(Z>=*IDA<2) 

M..1. J.P.l./.-M Mdiiflni an 



Jiisj-l 

LA=MA*JM 

LBbM8#JM 

MSU80190 

MSU80200 

MSU8C210 



LCMCMM 

MSUB0220 
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'Jp U 1 4 ». 



:WH 


?0 1 1 = 1 # ! A 
lAnL**i 

-IfakiAl— 

LCbLC^I 

C(LC)«A(tA)-B(LB> 

-R£ IURn 

End — — 

1 SUBROUTINE TO DECOMPOSE MATRIX FOR SIMULTANEOUS EQUATIONS 

CD60QDA + Q Q7 DATE 05 

decomTaTI 

DIMENSION A(MS1Z£,MSIZE), INTRCMSIZE) 

MATR IX DECOMP OSITION USED WITH SOL? SUBROUTINE FOR SOLUTION 

OF LINEAR SYSTEMS “ — ^ 

IF MATRIX A IS SINGULAR INTR(N) WILL BE SET TO .ZERO 

K = NN 

ITR=1 

NMsN-1 

DO ID J«1,NN 

AMAX«ABS<A( J, J) ) ” “ 

JPa J*l 


DO 2 IaJP,N ” 

AT-ABSCAC I# J) ) 


hi jrfil 


AMAXsAT 

INM 


I F < AHAXM,J,4 

3 lNTR(J)aJ 

GO TO 11 

IF ( IN)5,7,5 

MTRs-nTR 

DC 6 I s J_#_N 

A T*A ( J# J ) 

A<J, I )»A( IN, I) 

_*-U±A_Ll“ AT 

I N T R { J > a I N 
AnAXo-l.o/A(J,J) 
DO in 1= JP,N 
I F ( A ( I , J ) ) 8 , 1 0 » B 
8 AT=A( I, J)*AMAX 

AC I, J)*AT 

DO V KsjP,N 


I'iqdiLNll 
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9 

. lo 

M1,K>»A(J,K)#AT*A<I,K> 

CONTINUE 

— Lf IA(NiN)) 12, 11,12 __ 



DEC0H380 

DECOM390 


U 

12 

NTR*C 

I^TR(N)»KTR 

-JjgJURN 



DECOM410 

DEC0M420 

• BTMN 

ft 

END 

DETERMINANT EVALUATION 
■ — CD600D4.007 niTP 

SUBROUTINE 

09 >0 A /AM 

DEC0M44 0 
DTMN0020 

ft 


SUBROUTINE DTMN(A,1NTR,MSIZE*NN,DET) 
DIMENSION ACMS1ZE,MSIZE), INTRIMSIZE) 
GCmilES DET, THE DETERMINANT OF THE 

DECOMPOSED 

MiT&iy 

u T HJi U D J 0 ... 

0TM.NOO4O 

DTMN0050 

•ft 

ft 


subroutine decom must be called first 

INTR(N) WILL CONTAIN INTEGER POWER OF 

TEN OF MULTIPLYING FACTOR 

DTMN0070 

DTMNooen 



N E ® 3 8 

n*nn 

M*L»INTR(N> 

- 


OTMNOIOO 

DTMN0110 


1 

IF(NTR)2,1,3 
DTTsO.O 
-GO TO IQ 



OTHN0130 

DTMN0140 


2 

a_ 

DTTs-10, 

GO TO 4 

-U.Tafltl 



DTMNU160 

DTMN017P 


4 

00 9 1=1, N 
OTsAPS(DTT) 

_lf_LABS< At I^U )M .0 >9,0.7 



OTMN0190 

DTMN0200 


5 

6 

IFIDT-l. )6#9,9 

DTTaCTTftEP 

NTRsuTR-NE 



D TMN 0220 
0TMNO23O 

0TMNfl?4n 


7 

GO TO 9 i 

I r < DT-l . 0 )9,9,8 ; 

DTT.OTT/EP 



DTMN0250 

DTMN0260 


9 

NTR = NTR* N£ 
DTT»OTT*A( I, I ) 
iNTPf N)» nTR 



DTHN0280 
DTMN0290 
Dtmnci 3fi n 


10~ 

OET=0TT 

RETURN 

£ is li 



DTMN0310 
DTMN032 D 

*PHAT 

SUBROUTINE TO 

SUBROUTINE PMAT(A,NR,NC,MM,NN) 
DIMENSION A(MM,NN),P(6) 

PRINT MATRIX 

— R 1 n M U w •? U 

PMAT0020 



KROW=NR 



PMATQ060 


72 


£col«nc 
1*1 
J =1 

iTpm — 

JP = J 

_ D 0 _? K = 1 , 6 
KK-K 

p<k>*aci,j> 

— jaj* 1 

2 CONTINUE 

— CALL ne hun 

WRITE (6# 4 » I P # JP, <P<0,J(«1»KKT 
4 P0RHAT<2M,6(1PE16,8) ) 

QO TP 1 

3 CALL NENLIN 

*»°1 

IF< t.lE,NROW)QQ TO 1 


OETMOW/GETMOR — ~ — 

GETMOR 

GETMOR 

TO 08TAIN MORE CO^ToOTIc 

CALL GET.MomYPE, RESULT, HUM, FC> 

TYPEcfl f OR CORE 
TYPE*! TOR RANDOM DISC 
TYP6«2 fCR LINKED DISC 
-RESUL T* 0 IF SUCCESSFUL 

RESULTsi if UNSUCCESSFUL 

NUM IS NUM8ER OF LINKS DESIRED OR 

NUMBER OF K (1024 HOPiKf DE SIRED 

PC IS THE FILE CODE IN THE FORM 6H0QOOFC 
(USED ONLY WHEN GETTING MORE DISC) 

4*1* 

2,1# — 

* 6 H 0 | JO . OO 1 

- ST ORA CURE REQUEST 

R AND0M RANDOM DISC REQUEST 
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) 


r- " r " 



LDC 5,1* 

TR4 ORTOA 

RANDOM LAG § t ^» 

ORQ *1,DU 


LINKED DISC REQUEST 


ORTOA OR A =2#DU 


ZERO 


MME 

ZERO BSS 

m 

GEMORE 

1 

mo 



■6HOOQOQO 

CONT 

sfiHflQOOOl 


CONT STQ 

tra 

EAiD 

3*1* 

0*1 


i 

s 

ghap 
incode 
Ill 

DECK.COMDK 01524 Q90672GIMME 

IBMF 

—Gl-y.E-HE-HQRP-r.OPP fMMznnm 


L»L 

GIMME# GIVE-ME-MORE-CORE 

IMMEO0G2 

• 


J.Q ADJUST SIZE OF SPECIFIED ARRAY TO SPEC ! F 1 Fn 

IMM60003 

• 

• 

_ -t 

NfcW SIZE, MME GEMORE IS USED TO ADD MORE WORDS OF MEMORY, 
MME GEMREL IS USED TO RELEASE SURPLUS WORDS OF MEMORY. 

IMME0005 
I MME 0 0 0 6 

• 

• 

« 

CALLING SEQUENCE 

CALL GIMME ( NEWS1Z# OLDSIZ, ARRAY NAME ) 

IMMED008 
IMMEOOC9 
I MMFn n i n 

• 

• 

• 

THb ARRAY SPECIFIED MUST 0E LOADED AT THE VERY TOP OF MEMORY, 
USING EITHER A S USE CONTROL CARD OR A BLOCK DATA SUBPROGRAM, 
. I IT IS ALSO NECESSARY THAT THE ARRAY BE IN' LABELED COMMON. ) 

IMMEQOll 

IMME0012 

IMHE0D13 

• 

• 

• 

WARNING — THIS SUBROUTINE CHANGES THE SECOND ARGUMENT 
TO REFLECT THE NEW SIZE OF THE ARRAY 

IMME0014 
I MME0015 
I MMEO 01 6 

• 

SY*REF 
GIMME. SAVE 

.FCNV, 

IMMECQ17 
I MMEtiOla 

• 

• 

• 

these next 

THEY MAY BE 

INSTRUCTIONS ARE USED TO MONITOR PERFORMANCE OF GIMME 
REMOVED IF NOT WANTED# ALONG WITH. THE CODE AT SYMBOL 

I MMEO 019 
IMME0020 
1 MHE0021 

• 

• 

DOME. SET 



DONE EQU GIMME+1 IF REMOVED, 
-QE-TJHE . GET TIME OF DAY 

I MME 0 0 22 
I MME 00 23 


ST3 

TEMP2*1 

I MMEO 025 
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• 

mme 

STO 

GELAPS 

TEMP2 

GET PROCESSOR TIME USED TILL 

NOW 

IMME0026 

IMME0027 

IMMEQG28 

• 

UNfcUK THAT 

THE REFERENCED 

ARRAY IS LOADED AT THE TOP OF CORE 

IMME0O29 






fMM£G030 


SBAR 


GET THE ADDRESS OF THE TOP OF 

CORE 

IMME0031 


LDO 

*•* DL 



IHHEQ032 


ANO 

QLS 

■0777 / DL 
9 



IMME0033 


EAXO 


ffff ADDRESS OF AftftAV 




STXO 

ARADR 

SAVE 



ARADR S0O 

**» DL . 

SUBTRACT ADDRESS OF ARRAr 




LDA 

LOWLOD 

CHECK LOWLOAD 


a ... r If 


CAKA 

»1B10#DL 

INDICATOR 




tnz 

-LOW 

...LOWLOADED 



- * 
_ * 

MIGHLOAD PRE-PROCESS AND 

CHECKS 




LDA 

• 1 , D U 

SET FLAG 




STCA 

RLSFLO# 70 

TO ALLOW MEMORY RELEASE 




STO 

— T_E HP 

SAVE ARRAY SIZE 




SB Q 

3,1* 

SUBTRACT CLAIMED ARRAY SjZE 




TZ6 

*♦3 

...EQUAL 




CftPQ 

IjOU 

MAYBE ARRAY WAS LOADED On ODD 

LOCATION 

IMME0039 


Tf-iZ 

ERROR 

... SORRY, SOMETHING ELSE AT 

TOP OF CORE 

I M M c 0 0 4 (J 


LDO 

TEMP 

RESTORE ARRAY SIZE 



• 

TRA 

GIM 

...CONTINUE 




* LOWLOAD PREPROCESS AND CHECKS 


LOW LDXO ARADR " 'ADDRESS OF ARRAY 

CMPXO LIMITS COMPARE WITH LOWEST UNUSED 



tmi 

ERROR 

...ARRAY NOT LOADED ABOVE PROGRAM 



L X Hi 
CMPXO 

—III 

LIMITS 

ARADR 

• 44 

ADDRESS OF HIGH UNUSED LIMIT 
COMPARE WITH BASE OF ARRAY 
.i 4-PKj. LIMIT IS BELOW B A_Si 



LDXO 

S8LX0 

SXLC 

ARADR 
*1, DU 
LIMITS 

ADJUST LIMIT 
TO 8E 

BELOW ARRAY 



LDA 
SB A 
STCA 

2,1* 

3il* 

RLSFLG, 70 

COMPUTE FLAG 

INDICATING THAT 

RELEASE MAY OCCUR (IF NEGATIVE) 


gTR 

“5TC 

3,1* 

MAKE SURE OLE ARTAY SIZE IS CORRECT 
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T~ ■■r 



• 

• 

• 

COMPUTE 

THE number OF 1024 

WORD BLOCKS THAT ARE REQUIRED 

IMME0041 

IMME0O42 



LOO 
SBO 
- ADO 

2,1* 

3,1# 

— =1023,01 

NUMBER OF WORDS REQUIRED 

NUMBER OF WORDS ALREADY IN THE ARRAY 

ROUNDING UP FACTOR 

IMME0044 
• 1MME0045 



QRS 

STQ 

- T.ZJL- 

10 

TEMP 

D.Q.NE 

DIVIDE BY 1024 
SAVE NUMBER OF K TO GROW 
. ...NO CORF ADJUSTMENT REoUlRFn 

IHHE0050 

• 

• 

IF GIMME 
HERE 

is not to release 

CORE, INSERT THE FOLLOWING INSTRUCTION 



__ X?LL- 

DONE 

...NO MORE CORF REQIlfPEn 




LDA 
HS 
— LSA . 

TEMP 

10 

lxJL« 

UPDATE THE NUMBER 
OF WORDS IN 
THE ARRAY 


• 

• 

IF gimme is not to release 
INSTRUCTIONS. 

CMPA 1 LdLL 

SURPLUS CORE, REMOVE THE FOLLOWING TWO 

1MM6D056 

IMM6Q057 

• 


TH I 

GIVEUP 

... CORE CAN BE RELEASED 

IMHE0Q59 

• 

LOOP TO 1 

FETCH REQUIRED CORE 

..IN 2K INCREMENTS 

IMME0060 

• 

LOOP 

STQ 

CMPO 

TEMP3 

JUHL 

SAVE THE NUMBER OF K 
HOW MANY MORE K 

IMME0C62 




LAST 

GEMORE 

1x2 

... GET IK MORE 
OET 2K MORE 

I MMfc 0 0 6 4 
IMME0065 



TSX1 
LOO 
So Q 

CNT 
TEMP3 
2. DL 

... CORE REOUEST REFUSED 
UPDATE THE 
NUH8ER OF K 

IHME0067 

IMM60069 

LIST 

T2E 

tsa 

K F 

DONE 

LOOP 

GEMORE 

... don! 

... CONTINUE 

QET THE LAST IK BLOCK 

IMME007D 

IMME0072 

T MNPfl ft7Tt 



ZERO 

TSX1 

ISi^_ 

0,1 

CNT 

nn»j^ 

... CORE REQUEST REFUSED 
-x_. t DONE 

illlM 

• 

CNT 


AOS 

LOS 

COUNT 

■ 4»1S Q i*64 , DL 

COUNT NUMBER OF REQUESTS REFUSED 
GO TO SLEEP FOR A UMTIF 

IMME0077 
I M M E O O 7 8 

• 


MME 

tra 

GEWAXE 

-3,1 

... ZZ2 2 Z Z Z Z Z 

... TRY AGAIM 

IMMEQ080 

IMME0081 

* 

SEQUENCE 

TO RELEASE SURPLUS 

CORE 

IMME0083 
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-*• 

GJVEUP 

SZk 

**, DU 

CHECK RELEASE FLAG 

IMHE00B4 


TPL 

DONE 

...POSITIVE, NO RELEASE REQUESTED 



ICQ 

TEMP 

GET ABSOLUTE VALUE OF NUMBER OF BLOCKS 



QIS 

10 

MULTIPLY BY 1024 

1MME0086 


EAA 

DONE 

SET RETURN ADDRESS FOR GEMREL 

JMME00B7 


M M E 

GEMREL 

... RELEASE CORE 

IMME0088 

* PRINT STATISTICS FDR THIS 

CALL TO GIMME 

I MMEC089 
IMME0093 





IMME0091 

DONE 

MMF 

6ELAPS 

COMPUTE PROCESSOR TIME USED BY GIMME 

IMME0092 


SSQ 

TEMP2 


IMME0093 


MME 

GETIME 

COMPUTE ELAPSED TIME 

INME0094 


SPQ 

TEHP2*1 


IMME0095 


TPL 

• + 2 

...DID NOT PASS MIDNIGHT 



ADQ 

=5,5296E9 

ADD 24 HOURS TO CORRECT ' 



TOV 

**1 

CONVERT TO FLOATING POINT 

IMME0096 


LDA 

0,DL 


IMME0097 


IDE 

■71325, DU 


IMHE0098 


FNO 



IMM60099 


FD V 

■2.304F8 

CONVERT TO HOURS FROM 64THS OF MILLISEC. 

IMMF0100 


FST 

T6MP2*1 

SAVE FOR PRINTING 

IMMtOim 


ldq 

TEh°2 

CONVERT PROC. TIME TO HOURS 



TOV 

*♦1 




LDA 

0 , DL 




LDE 

*71025, DU 




FNO 





fdv 

■ 2.3ME8 




FST 

TEHP2 




CALL 

.FPRN»(FILE,FORMl) 



LDA 

TEMP 

PRINT NUMBER K CORE 

1 


TSX1 

• FCri V » 




LDA 

COUNT 

_ .PRINT NUMBER OF REOUESTS REFUSED 

I MME0103 


TSX1 

•fcnv* 


I MMEO 1 05 


LDA 

TEMP2 

PRINT PROCESSOR TIME USED BY GIMME 

IMME0106 


TSX1 



IMMEQ107 


LDA 

TEMP2*1 

PRINT ELAPSED TIME 

immeoiub 


TS XI 

.fcnv. 


IMME01D9 


CALL 

.FFIL. 


IMMEQllO 


STZ 

COUrT 

RESET FOR NEXT CALL TO Q I HME 

I MM£0 111 


RETURN 

GIMME 


IMMEQ112 





IMME0113 

ERROR 

C4LL 

FXEM ( *61, MESS, 

■ 5) " " 



77 



rlsflg 

ECU 

GIVEUP 



LIMITS 

bool 

37 



LOKLO! 

bool 

—2A 



TEMP 

8SS 

1 


I MME0115 

TEMP2 

BSS 

2 

STORAGE FOR PROCESSOR AND ELAPSED TIME 

IMME0U6 

1 

B SS 

-1 



COUNT 

ZERO 



I MME0117 

MESS 

BC I 

5, ARRAY NOT LOADED ABOVE PROGRAM 







3CI 

8, USED # F7*6# 11H HR, PRQC. # # F7, 4# 11H HR. ELAPS.) 



block 

GIMMES 



— FILE . 

DEC 

06 




END 



IMME0121 

S 

OMAP 

DECK# COKDK 

77739 C11671T I ME 

s 

INCOPE 

IfiMF 




LBL 

TIME, TIME 

AND DATE SUBPROGRAM 

T1ME0010 


TTL 

DRILL CONVERSION PROGRAM 

TIME 0020 


-HU— 

TIME* DATE 



• 




T I ME 0 0 40 

• CALL TIMDAT (TINE f DATE> 

• 

TIMEOO50 
T t mpo n*n 

WHERE TIKE 

= 2 CONSECUTIVE WORDS WHERE THE CURRENT TIME WILL RE 

TIME0070 



PLACED AS 

HHJMMJ SS 

TIME008O 

« 

DATE 

= 2 CONSECUTIVE WQRnS WHERE THE CURRENT DlTP WtIL fif 




PLACED AS 

MM/DD/YY 

T I MEOl 0 0 

* 




TIMEQ110 

TIMDAT 

SAVE 

0 


TIME0120 


NHE 

GETIME 

GET DATE IN A AS MMDDTY AND TIME IN Q IN 

T I ME 0130 


REM 


64*THS MSEC SINCE MIDNIGHT, 

T I ME014 0 


ST Q 

TIME 

SAVE TIME 

TIME0150 


LRL 

36 

PLACE DATE IN Q 

T IME0160 


LDA 

*6H 

PLACE SPACES IN A 

TIME0170 


llb 

_ X2 

MOVE MM INTO A 

T I MEOlfl n 


als 

6 


TIME019D 


ORA 

=3hgo/,dl 

INSERT SLASH 

TJMEQ2Q0 


LLR 

i-2 _ 

MOVE DD INTO A 

Tj ME 0210 


ALS 

6 


T I ME 0 22 0 


ORA 

33HOO/#DL 

INSERT SLASH 

TIME0230 


LUXC 

Itl 

LOAD ADDRESS OF DATE 

T I ME 0 2 4 0 


STA 

0,0 

STORE FIRST 

T I ME025 0 


STQ 

1,0 

AND SECOND WORD 

T I M E 0 2 6 0 



time ■ 

PUT TIME IN 0 

T IME0270 


ORS 

6 

CONVERT TO MSEC 

TIME0280 
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D I V 
DIV 
STA 

1000 # DL 
10, DL 
time 

CONVERT TO SEC 
STORE UNIT SECONDS 

TIHE0290 

TJME0300 

TIME0310 


DIV 

STA 

DIV 

6,DL 
TIME+1 
10, DL 

STORE TEN’S OF SECONDS 

TIME0320 
T I ME0330 
TIME0340 


STA 

DIV 

STA 

time*z 

6# DL ‘ 
T1ME*J 

STORE UNIT MINUTES 
STORE TEN’S OF MINUTES 

TIME0390 
• TIME0360 
TJHED370 


DJV 

ALS 

LLR 

“lo# D L 
30 
6 

UNIT HRS TO' A,. ID’S OF HRS IN fl 

ASSEMBLE HH IN RIQHT OF Q 

TIME0380 
T I MEQ390 
Tl ME 040 0 


QLS 

ORO 

6 

*015 , DL 

INSERT COLON 

T I ME0410 
T I ME 0420 


- qls 

6 


TrMEfUSfl 


ORQ 

QLS 

ORO 

TIMERS 

6 

T I ME*2 

INSERT MM 



QLS 
ORQ 
— LDA 

6 

=015, DL 
= 6H 

INSERT COLON 
PLACE SPACES IN A 

TIME0470 

TIME04B0 

T I ME Q 49 n 


LLS 

ORQ 

—LLS 

6 

TIME*1 

(l 


T I ME 0 5 0 0 
TIME0510 
T T M£052n 


ORO 

LLR 

Lnxo 

TIME 

24 

JU1 

INSERT SS 

LOAD ADDRESS OF T T ME 

T I M60530 
TIM60540 
Tf MEHS«5rt 


STA 

STQ 

RETURN 

0,0 

1,0 

T IMP AT 

STORE FIRST AND 
SECOND WORD 

TIME0560 

TIME0570 

TIME0580 

• 

• 

EJECT 

CALL ELT I ME { T I MEL ) 


TIME0590 

TIME0600 

TlMEOfilO 

• 

* WHERE TIMEL 

* LOCATION 

where elapsed time in msec, is placed. 

TIME0620 
time 0630 





TIME0640 

ELTIME 

SAVE 

MME 

QRS 

gelaps 

6 

set ELAPSED TIME IN Q IN 64*THS MSEC. 
CONVERT TO MSEC 

TIME0650 

TIM60660 

TIME0670 

TIME 

STQ 

RETURN 

BSS 

2,1* 

ELTIME 

4 

STORE IN TIMEL 

TIME0680 
T I HE 0690 
T I ME 070 fl 


EXECUTE 
FFILE (LI 


S 

I 


FILE 08,,1R 
LIMITS 2QQ,17K,,20X 





